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Foreword 


This report presents the results of a literature review as part of the development of composite 
hardware fracture control guidelines funded by NASA Engineering and Safety Center (NESC) 
under contract NNL04AA09B. J. B. Chang is the principal investigator. 

The objectives of the overall development tasks are to provide a broad information and database 
to the designers, analysts, and testing personnel who are engaged in space flight hardware 
production. 

The literature review concentrated on the state-of-the-art damage tolerance analysis methods that 
have the potential to be used in the damage tolerance demonstration for fracture-critical 
composite hardware in manned and unmanned spaceflight systems. 
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1.0 Introduction 


The National Aeronautics and Space Administration (NASA) requires that all payload that used 
the Space Shuttle and all equipment installed in the International Space Station comply with the 
fracture control requirements specified in NASA-STD-5003, Fracture Control Requirements for 
Payloads Using the Space Shuttle and SSP 30558B, Fracture Control Requirements for Space 
Station. In a general fracture control standard, NASA-STD-5007, General Fracture Control 
Requirements for Manned Spaceflight Systems, fracture critical items (FCIs) shall be 
demonstrated to have adequate damage tolerance capabilities, i.e., to have adequate “safe-life.” 
This requirement is applicable to all FCIs independent of materials. 

NASA-STD-5007 has been replaced by NASA-STD-5019, Fracture Control Requirements for 
Spaceflight Hardware. In this standard, fracture control of the composite and bonded structures 
specified in MSFC-RQMT-3479, Fracture Control Requirements for Composite and Bonded 
Vehicle and Payload Structures, shall be applied to all fracture critical composite hardware. The 
detail requirements specified in these standards are shown in NASA-STD-5019. 

In fracture control for hardware items made of metal, the primary failure mode that should be 
considered is the crack propagation leading to catastrophic failure. Cracks or crack-like defects 
are the “damage” to be assessed. Damage tolerance demonstration consists of life and residual 
strength predictions that are obtained by performing fatigue crack growth analysis and static 
fracture strength analysis. Linear elastic fracture mechanics (LEFM) principle is the base for 
such analyses. The primary parameter used in LEFM is the crack-tip stress intensity factor, 
which is a function of the crack size, crack geometry and applied loading conditions. Crack 
growth analysis computer codes used in the aerospace industry, including CRKGRO (Chang 
et al., 1981) and FLAGRO (Hu, 1980), all contain various Absolutions and a large database of 
fracture properties (K c ) and crack growth rates (da/dN). The capabilities of crack-growth life 
predictions using LEFM have been shown to be adequate when correlated with experimental 
results (Chang, 1981). NASGRO (JSC-22267B, Fatigue Crack Growth Computer Program 
NASGRO®, 2001), a derivative code from FLAGRO, is widely used for space flight hardware. 

In fracture control of composite hardware, it is generally agreed that analytical methods that can 
be used for damage tolerance demonstration of composite structures are not as mature as the 
metallic counterparts. It is known that tensile failures of composite structures are primarily K 
associated with fiber breaking. On the other hand, compression and shear failures of composite 
structures are often associated with matrix cracking and local delamination. Furthermore, a 
commonly used composite material system such as graphite/epoxy (Gr/Ep) is very brittle and 
highly susceptible to impact damage. Therefore, the evaluations of damage tolerance capabilities 
of composite hardware are much more complicated than their metallic counterparts. 

To meet fracture control requirements levied in NASA-STD-5003 and SSP 30558B, it is 
common practice in the spacecraft industry to perform acceptance proof tests on all flight 
composite structures at a 1.2 x limit load level. For many applications, it has been shown to be a 
viable option to meet the requirement set forth in those standards. However, unless there is 
reasonable compensation for the effect of environment, proof tests are required to be conducted 
at the service environment. In some cases, to perform a proof test at a specific environment is 
both costly and time consuming. For example, a composite cryogenic tank adapter used in an 
upper stage vehicle should be proof tested at cryogenic temperature. However, due to facility 



limitations and other technical concerns, it may be decided to conduct a damage tolerance test 
program instead (Sollars, 1987). Furthermore, to avoid a high proof load causing matrix cracking 
and fiber breakage, there is a limit for a proof test level for composite structural items. NASA- 
STD-5003 specifies that the proof test level should be limited to less than 80 percent of ultimate 
strength. For a marginal design, the proof load should be <1.1 x limit load; in this case, the 
“1.2 x limit load” requirement cannot be met. Hence, damage tolerance testing is the only viable 
option if damage tolerance analysis is not acceptable by most fracture-control authorities. 

The aircraft industry has traditionally used the “building block” approach to demonstrate the 
damage tolerance capability of fracture critical airframe structures (Whitehead and Dee, 1983). 

In the building block approach, one would run different level tests in concert with damage 
tolerance analysis. This approach is currently being adopted by the space industry for the same 
purpose (Konno et al., 2001). Because there are many impact threat events, such as tool drop, 
runway debris, hail, lighting strikes, etc. in the operation of an aircraft, impact damage is the 
predominant “defect” to be considered. Compressive strength after impact is used as the 
allowable for composite wing and other critical airframe structural designs. Many attempts have 
been made by investigators to predict residual strength for those structures (Rhoades et al., 1978; 
Byers, 1980; Smith and Wilson, 1985; Horton et al., 1988; Madan, 1988). Another important 
defect that needs to be evaluated is delamination. The growth of delamination in compression- 
compression fatigue spectrum applied structures such as upper wing covers was the focus of 
research in the 1970s and 1980s (Rybicki and Kanninen, 1977; Chai and Babcock, 1985). From 
this research, a few growth models have been proposed. 

For spaceflight composite structures such as satellite trusses and platforms and launch vehicle 
fairings, de laminations are very critical defects since these structures are also subjected to 
compression loads. For other composite hardware items such as overwrapped pressure vessels 
and pressurized structures, impact damage and other types of mechanical damage including 
surface cuts and scratches are potential damages that need to be considered. Damage tolerance 
testing on some of the above-mentioned composite hardware could be very costly. Thus, the use 
of damage-tolerance analysis technologies to demonstrate the damage tolerance ability of 
composite FCIs should be encouraged. Accordingly, this state-of-the-art review was conducted; 
the results are presented in this document. 

2.0 Residual Strength Prediction Methodologies 

Damage tolerance analysis methods reviewed in this task cover two categories: residual strength 
analysis and “safe-life” analysis. Defects included are surface cuts and/or scratches, delamination 
and impact damage. The residual strength analysis methods are presented in this section while 
the life prediction methods are discussed in Section 3. 

2.1 Surface Cuts and Scratches 

2.1.1 Closed Form Solutions for Surface Cuts and Scratches 

For laminated composite structures, it has been observed by many investigators that defects such 
as cuts, scratches and impact damage will in general not grow under tensile cyclic loading. In 
these cases, damage tolerance assessments require the consideration of the residual strength 
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instead of “safe-life.” A few methods and models have been developed for predicting the 
residual strength of these types of damage. 

• LEFM (Classic) Method 

• Inherent Flaw Model 

• Point Stress Model 

• Point Strain Model 

• Mar-Lin (ML) Model 

2. 1.1.1 Brief Description of Closed Form Solutions 

LEFM (Classic) Method 

The LEFM method uses the crack-tip stress intensity factor, K, to characterize the crack behavior 
in an isotropic body. This method has been routinely used to predict the residual strength of 
cracked metallic spaceflight structures since the 1970s. For a through-the-thickness cut in a 
quasi-isotropic composite laminate plate subjected to tensile load, attempts have been made to 
use the same parameter, K. In this case, the cut is treated as a crack. 

Based on the LEFM method, the strain in the fiber direction at a distance r directly ahead of a cut 
(crack) tip can be written in the following infinite series (Poe, 1983). 

£- 1 =Q(2^r)”X+0(r°) (2.1) 


where 


Q=* c <r/E x 


( 2 . 2 ) 


and 


£ = (i - V ) ( V E x/ E y sin2 ^ + c °s 2 «) ( 2 - 3 ) 

In Equations (2.2) and (2.3), K c is the critical stress intensity factor (or fracture toughness) of the 
composite laminate, x and y are Cartesian coordinates with x perpendicular to the cut; E is 
modulus of elasticity, vis Poisson’s ratio, and eris the angle that the fiber makes to the x-axis 
(perpendicular to the crack). 

Inherent Flaw Model 

The inherent flaw model (IFM) was originally developed by Waddoups et al. (1971) (Aronsson, 
1986) for predicting the residual strength of notched laminates with infinite width under tensile 
load. It has been extended to finite width by introducing a finite- width correction factor. For a 
composite laminate with a notch (crack) of length 2a and inherent damage zone size, c 0 , 
symmetrically on both sides of the crack tip, the notched tensile strength cr™ is given by 

= (2.4) 

Y V a + c 0 
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where <J 0 is the unnotched tensile strength and Y is a finite-width correction factor. 

While the material is considered homogeneous, orthotropic and linear elastic, the following 
relation-ship between <j^ FM and the critical strain-energy release rate (fracture toughness), G c , is 

obtained (Aronsson, 1986) from the equation: 


_IFM 

°N 


Y Jl + ^- 

V CG„ 


(2.5) 
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66 


11 


( 2 . 6 ) 


In Equation (2.6), Ay are the in-plane laminate stiffnesses as determined from laminate plate 
theory (Ashton and Whitney, 1970). A linear relationship between c 0 and G c may further be 
obtained from 


c 


O 



(2.7) 


Point Stress Model 

The Point Stress Model (PSM) was derived by Whitney and Nuismer (WN) who considered the 
exact elasticity solution of the normal stress ahead of a notch. By approximating the exact 
solution of the normal stress ahead of the notch (crack) with the asymptotic solution, the model 
can be extended to notched laminates with finite geometry. This is justified if the notch length is 
sufficiently large compared to the size of the damage zone size d 0 (Agarwal and Giare, 1982; 
Zweben, 1973). 

For a tension loaded laminate with a notch length 2a, an approximate expression for the 
unnotched tensile strength according to the Whitney and Nuismer Point Stress Model (WN- 
PSM) is given by (Aronsson, 1986): 


_psc 



( 2 . 8 ) 


where do is the characteristic length defined in the PSM (Whitney and Nuismer, 1974; Nuismer 
and Whitney, 1975). 

The finite-width correction factor, Y, used in the calculations was obtained from (Brown, 1970): 


Y 



1.09-1.73 




(2.9) 
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This factor was derived originally for an isotropic material and is used for the quasi-isotropic 
laminates considered here. 


The WN-PSM (Nuismer and Whitney, 1975; Whitney and Nuismer, 1974), predicts failure when 
the stress at a characteristic dimension, di, ahead of the crack tip equals or exceeds < 7 0 . The 
residual notched strength is given by 


1 


1 - 


a + d 


1 J 


( 2 . 10 ) 


The two parameters in this model, <J 0 and di, are usually determined by tests. 

PSM 

The PSM proposed by Poe and Sova (Poe, 1983; Poe and Sova, 1980), may also be formulated 
with a characteristic dimension, d 2 . This model predicts failure when the strain at a distance 
ahead of the notch (crack) tip equals or exceeds the fiber failure strain. The notched failure stress 
is given by 



( 2 . 11 ) 


where £ is a functional that depends on elastic constants and the orientation of the principal load 
carrying plies. The characteristic dimension relates to a material toughness parameter, which was 
found to be relatively independent of lay-up. The two parameters that must be determined for 
this model are the fiber failure strain and d 2 . 

ML Model 

The Mar-Lin (ML) model (Lin and Mar, 1977; Mar and Lin, 1977) allows the singularity, n, in 
the stress intensity factor to be a value other than square-root. The notched failure stress is given 
by 


= 


H„ 


(2a) n 


( 2 . 12 ) 


where H c is the composite fracture toughness. In general, H c and the exponent n are the two 
parameters that must be determined by test. In studies, the exponent, n, was related to the 
theoretical singularity of a crack in the matrix, with the tip at the fiber/matrix interface. In this 
case, the singularity is a function of the ratio of fiber and matrix shear moduli and Poisson’s 
ratios. Using this method, the singularities for a range of typical fiber/matrix combinations were 
determined to be between 0.25 and 0.35. 

2. 1.1. 2 Comparison of Various Closed Form Models 

The following discussions are taken from CMH-17, Composite Materials Handbook Volume 3. 
In this reference, predictions for both small notch (2a ~1.2 in. (3 cm)) and large notch (2a up to 
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20 in. (0.5 m)) sizes are compared. The former notch sizes are characteristic of much of the data 
collected for composites to date. Four models are compared; LEFM (Classical), WN (point 
stress), point strain (PS), and ML. As a baseline for comparing changes in notch length predicted 
by the four models, curves are generated based on average experimental results (finite width 
corrected) for the IM6/937A tape material with W/2a = 4 and a lay-up of [+45/90/-45/0/+30/- 
30/0/ — 45/90/+45]. This will ensure that all theories agree for at least one crack length. 

Figure 1 shows a comparison of the four models for small crack notches. Only a small difference 
is seen between PS and WN models. A close examination of the LEFM and ML curves indicates 
that the singularity has a significant effect on curve shape. For notch lengths less than the 
baseline point, ML predictions are less than those of Classical. For notch lengths greater than the 
baseline point, the opposite is true, and models tend to segregate based on singularity, i.e., WN, 
PS and Classical yield nearly the same predictions. 

Figure 2 shows that singularity dramatically affects differences between predictions in the large 
notch length range. The ratio of notched strength predictions for models with the same order of 
singularity becomes a constant. For example, WN and LEFM become functionally equivalent, 
and the relationship 

K lc =a 0 ^d[ (2.13) 

will yield a value for K\c such that the two models compare exactly for large notches. 



Figure 1. Comparison of curve shapes for notched strength prediction theories in small notch range. 
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Figure 2. Comparison of curve shapes for notched strength prediction theories in large crack notch. 
2.1.2 Finite Element Analysis (FEA) Methods for Surface Cuts and Scratches 
The FEA method provides flexibility and accuracy for multiple configurations. 

2. 1.2.1 Conventional FEA Methods 

Two conventional FEA methods exist to account for the effects of damage progression on load 
redistribution in finite element models (FEM). Progressive damage methods that degrade various 
stiffness properties of individual elements as specified failure criteria are met (Dopker et al . , 
1995) have shown some success in modeling damage growth in specimen configurations. The 
magnitude of the calculations, however, provides a significant obstacle to incorporating them 
into the complex models required for stiffened structures. 

Notched Strength Modeling Using FEA 

In order to achieve a truly general notched strength model, the physical basis of failure must be 
considered, and an extensive study of subcritical damage is an essential element of this process. 
To this end, a notched strength model was developed by Kortschot and Beaumont (1991) based 
on experimental results that illustrate the dependence of notched strength on subcritical damage. 
The model was finite element based and used to determine the modified notch tip stress 
distributions. The failure criterion was complemented by theoretical predictions of damage 
growth. 

The approach of this method is to calculate the stress distribution near the notch tip for a FEM 
that includes a subcritical damage region. Using stress results from the FEA, the stress 
concentration factor, K t , can be determined. The failure criterion is derived from curve fitting the 
finite element results to test data. 

An example of this method uses a one-quarter model of the specimen using appropriate boundary 
conditions, as shown in Figure 3. Two-dimensional (2-D), eight-noded quadrilateral plane stress 
elements were used. Two layers of elements were superimposed with one layer of elements 
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given the properties of a 0° ply and the other given the properties of a 90° ply. Corresponding 
elements in the two layers shared all nodes and hence deformed identically except in the area 
representing the delamination region where the 90° ply was disconnected from the 0° ply. In 
addition, the elements of the 90° ply were continuous across the split in the 0° ply, as illustrated 
in Figure 4. Although this geometry is quite complicated, it is a direct representation of the 
observed subcritical damage in [90/0] s laminates. 

The FEM was constructed to yield information about the tensile stress distribution in the 0° ply 
near the notch tip. The model produced stress contour maps that displayed finite maximum 

stresses for all values of Ha 0, where t is the split length and a is the crack length (V-notch). 
The stress concentration factor, K t , was determined as the maximum stress in the 0° ply divided 
by the remote stress on the laminate for each Ha ratio. The relationship between K t and Ha was 
approximated very closely by the empirical function 

K t = 8.16(f/a)"°' 284 (2.14) 



□ 0° ply 
1. J 90° ply 


Figure 3. Finite element mesh for f/a = 1. 
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Figure 4. Schematic of typical subcritical damage in a [90/0] s laminate. 

In order to derive a failure criterion, that is, to discover what aspect of the 0° ply stress field is 
constant at failure, the terminal stress distributions for a variety of specimens must be compared. 
To do this in a simplified way, the peak stress in the 0° ply immediately prior to failure, Oq p , was 
calculated for all specimens. A value for terminal f/a was obtained from a radiograph of the 

failed specimen and substituted into Equation (2.14). The remote stress at failure, o^f, was then 
multiplied by K t to obtain Oo p . The predicted failure stress can be expressed as 

&oof( predicted) = (2.15) 

K t 


where 

K t = 8.16(f/a)-0.284, 

<7oof = remote failure stress of the laminate; and 
Oof = failure stress of 0° ply. 

Strain Softening Models 

Strain softening models (Llorca and Elices, 1992; Mazars and Bazant, 1988) have been 
successfully used to simulate the fracture of concrete and other materials by providing the ability 
to capture the global load redistribution that occurs as the notch-tip region is softened by damage 
formation, without the computational concerns of detailed progressive damage models. A range 
of softening laws has been proposed. However, determination of the strain softening material 
laws through trial and error requires a large number of tests, and can be computationally 
intensive. Basham et al. have presented an approach to determine these laws using energy 
methods and relatively few tests (Basham et al., 1993). 

This method is a generalized continuum approach that is compatible with complex FEMs 
required to properly approximate structural configurations. The approach allows the capture of 
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load redistribution due to local damage formation and growth. One approach is the use of 
nonlinear springs that undergo reductions in bending stiffness as damage occurs. The models can 
be calibrated using small-notch test results, and then extended to large-notch configurations. 
Modeling and calibrating the stiffness reductions are of concern for most structural 
configurations, where out-of-plane loading, load eccentricities, and bending loads are common 
(Kongshavn and Poursartip, 1999). 

Non-local Damage Theory Using Strain Softening 

Nahan and Kennedy have used the damage induced strain softening approach for degradation of 
laminate properties to simulate crack growth in laminate plates (Nahan and Kennedy, 2000). 

This approach defines bulk plate damage as a linear variation of orthotropic damage through the 
plate thickness. A damage term is defined that can range from zero, being a no damage 
condition, to one, which corresponds to complete material failure. A bulk laminate stress-strain 
softening relationship is derived and incorporated into a FEM of the laminate. The FEA will 
generate a simulated load history, in which the damage can be observed to develop at the notch- 
tip and propagate to eventual failure. The theory performs well in predicting tension fracture 
over a wide range of notch sizes as shown in Figure 5. 



Figure 5. Non-local damage theory predictions versus tested strength. 

Predicting Crack Growth Using Micropolar Strain Softening 

In another study, Kennedy used micropolar strain softening for predicting damage growth in 
composite plates with different notch geometries (Kennedy, 1999). The micropolar approach 
differs from classical elasticity theory in that the deformation of the body is described by a 
displacement vector and an independent microrotation vector (Eringen, 1968). The theory 
prevents singularities due to localization of deformation in strain softening materials under 
simple stress states. 

An eight-node, isoparametric quadrilateral element, as shown in Figure 6, was developed for 
carrying out the plain strain micropolar elasticity analysis. The element is similar to that used in 
classical elasticity theory, except the addition of a rotation degree of freedom, at each node. 

The theory was used to predict damage growth in a rectangular laminate plate with notches 
consisting of a circular hole, an elliptical hole, and a sharp crack. In all cases, the results showed 
that the specimen size had a relatively weak effect on the failure load predicted, which indicates 
that the usefulness of this model is limited to materials with this characteristic. 
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Figure 6. 2-D eight-noded isoparametric quadrilateral element. 

Finite Element Solution Issues 

Finite element solutions for problems involving strain softening material laws involve a number 
of complexities not often associated with static structural analysis. Specifically, singular stiffness 
matrices are encountered when material failure is occurring in a sufficiently large area and/or 
when the structure is buckling. A solver with a variety of robust, nonlinear solution algorithms 
that is capable of modeling strain-softening response for orthotropic materials should be selected. 
Arc-length methods, such as Riks (1979), have proven useful in dealing with snap-through 
stability problems, and were useful in initial efforts to model tension loaded strain softening 
problems. 

Solving the problem dynamically minimizes a number of numerical difficulties because damping 
and inertial forces can smooth out system response and greatly reduce numerical noise in the 
solution process. As the maximum load is reached, local failure occurs, thus accelerating parts of 
the system, and the numerical integration in time can be stopped when a minimum acceleration 
related to system failure has been achieved. This has proven to be a very accurate failure 
criterion for compression- loaded structural systems (Dopker et al., 1995). 

Element Size and Formulation 

Strain-softening laws and the finite element size are interrelated, due to the effect of element size 
on notch-tip strain distribution. A great amount of mesh sensitivity has been found when strain is 
incorporated into FEAs based on classical plasticity theory. The sensitivity occurs because as the 
mesh is refined, there is a tendency for the damage zone to localize to a zero volume, which can 
lead to the physically impossible condition of structural failure with zero energy dissipation. It 
has been found that larger elements (i.e., > 0.20 in. (0.5 cm)) result in less-severe, but broader, 
stress concentrations and in strain-softening curves with steeper unloading segments (Scholz 
et al., 1997). Element size, therefore, is another parameter in the strain softening approach that 
must be determined. Fortunately, damage in composite materials typically localizes on a 
relatively large scale, e.g. relative to plastic yielding at a crack tip in metal. Element sizes 
required to accurately predict notch-length and finite-width effects in compression are typically 
larger than those required for tension. 

Element formulation and strain softening laws are also interrelated. Limited studies of 4-, 8-, and 
9-noded shell elements found that higher order elements lead to higher fracture strengths and 
large damage zones (Dopker et al., 1994). 
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2.2 Delaminations 

Delamination, often referred to as interlaminar fracture, is one of the few instances where 
fracture mechanics formalism is applicable to fiber-reinforced composites on a global scale 
(Anderson, 1991). A delamination surface can be treated as a crack, where the resistance of the 
material to the propagation of the delamination is the critical strain-energy release rate (fracture 
toughness), G c . Since the delamination typically is confined to the matrix material, e.g., between 
plies, continuum theory is applicable, and the delamination growth is self similar, i.e., keeping 
the same aspect ratio. For composite structures under compressive loads, one of the most 
predominant failure modes is the growth of delamination. This type of “crack” has been studied 
extensively by the aircraft industry since the aircraft wing and other control surfaces are 
generally loaded in both tension and compression. 

2.2.1 Closed Form Solutions for Delaminations 

Closed form solutions are often used in the determination of the critical strain-energy release 
rates, G c , of a composite laminate specimen. Typical specimens include the double cantilever 
beam (DCB) specimen, the end-notched flexure (ENF) specimen and the edge delamination 
(ED) specimen. An advantage of the DCB specimen configuration is that it permits 
measurements of Mode I, Mode II, or mixed-mode fracture toughness (see Figure 7). The ENF 
specimen has essentially the same geometry as the DCB specimen, but the latter is loaded in 
three-point bending, which imposes Mode II displacements of the crack faces. The ED specimen 
simulates the conditions in an actual structure - tensile stresses normal to the ply are highest at 
the free edge; thus, delamination zones often initiate at the edges of a panel. 

For pure Mode I loading, elastic beam theory leads to the following expression for energy release 
rate (ERR): 


G i 


B E I 


where B is the beam’s thickness, and 


E 


t _2P[a 3 

3-Aj 


(2.16) 


(2.17) 
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Figure 7. Mode I, II and mixed-mode loading of DCB specimens. 
Assuming linear beam theory, the corresponding relationship for Mode II is given by 

r 3 P n a2 

11 4 B E I 


(2.18) 


Mixed-mode loading conditions can be achieved by unequal tensile loading of the upper and 
lower portions of the specimens. The applied loads can be resolved into Mode I and Mode II 
components as follows: 

Pl=|P L | (2.19a) 


Pu- 



ll. 19b) 


where Pu and Pl are the upper and lower loads, respectively. The components of G can be 
computed by inserting Pi and Pu into Equations (2.16), (2.17) and (2.18) above. 
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Linear beam theory may result in erroneous estimates of G, particularly when the specimen 
displacements are large. In this case, the area method provides an alternative measure of G. 
Figure 8 schematically illustrates a typical load-displacement curve, where the specimen is 
periodically unloaded. The loading portion of the curve is typically nonlinear, but the unloading 
curve is usually linear and passes through the origin. G can be estimated from the incremental 
area inside the load-displacement curve, divided by the change in crack area: 


B Aa 


The Mode I and Mode II components of G can be inferred from the Pi - Ai and Pn - An curves, 
respectively. 

(Chai, 1982). The one-dimensional (1-D) thin film model is shown below. 



Figure 8. Schematic load-displacement curve for a delamination toughness measurement. 

2. 2. 1.1 Thin-Film Models 

Three delamination models derived by Chai in the early 1980s include the 1-D Thin-Film Model, 
Thick-Colu mn Model and the 2-D Elliptical Model (Chai, 1982). The 1-D thin-film model is 
shown below. 

Consider the three stages in thin-film delamination and buckling as shown in Figure 9. The 
delaminated film of thickness h and length £ is part of an infinitely thick body characterized by 

Young’s Modulus, E, and Poisson’s ratio, v. The strain-energy release rate, G, can be expressed 
as 


Eh(l-v 2 ) 

G = f(£ 0 _£ cr )(£ 0 _ 3 £ cr ) (2.21) 

where So is the compressed strain and s cr is the critical buckling strain. 


14 




Shortening € x = ” € o 



The thin- film model provides a simple means of predicting delamination growth, but it has very 
limited application due to its simplifying assumptions such as infinite thickness and through the 
thickness delamination. 

2.2.2 FEA Methods for Delaminations 

The FEA method provides flexibility and accuracy for multiple configurations. Analyses 
methods have been presented that utilize conventional continuum elements. Other methods use 
singular elements that contains singularity to predict delaminations. 

2.2.2. 1 Conventional FEA for Delaminations 

Analyses that use conventional continuum elements, such as, quadrilateral and triangular shell 
elements or hexagonal and pentagonal solid elements, are referred to as conventional methods. 
One such method for predicting delaminations in a composite structure using conventional 
elements is the Virtual Crack Closure Technique (VCCT), first presented by Rybicki and 
Kanninen in a classic paper (Rybicki and Kanninen, 1977). 

Calculation of Stress Intensity Factor by FEA 

In their paper in 1977, Rybicki and Kanninen presented a method for evaluating stress intensity 
factors using a modified crack closure integral and a constant strain finite element analysis. 

The method uses Irwin’s crack closure integral. The assumption is that if a crack extends by a 
small amount Ac, the energy absorbed in the process is equal to the work required to close the 
crack to its original length (Irwin, 1958). The following is the mathematical interpretation 



where G is the energy release rate, cr y and r xy are the stresses near the crack tip, u and v are the 
relative sliding and opening displacements between points on the crack faces, and Ac is the crack 
extension at the crack tip. The above equation can also be expressed as G = G\ + G\\. 

Taking the limit and using the stress intensity factors K\ and K\\ gives 

2 j^2 

G l = and G u =-V-J3 (2.23-2.24) 

h h 

where y 3= 1 for plane stress and 1 - v 2 for plane strain and E is Young’s modulus. 

Rybicki used constant strain elements near the crack tip, and found a good match comparing 
with the other results, using more complex methods. The finite element approach for evaluating 
Ci and Gw from Equation (2.22) is based on the nodal forces and displacements as shown in 
Figure 10. 

In equation form, the work required to close the crack an amount Ac is 

G i= lim T^-Fc-(v c -v d ) and G n = lim -2— T c (u c -u d ) (2.25-2.26) 

Ac->0 2Ac Ac^O zAc 

where the value of F c and T c are the y- and x- forces, respectively, that are required to hold nodes 
c and d together. 
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Figure 10. Finite element nodes near crack tip. 


The energy domain integral methodology (Shih et al., 1986; Moran and Shih, 1987) is a general 
framework for numerical analysis of the /-integral, energy release rate. This approach is 
extremely versatile, as it can be applied to both quasistatic and dynamic problems with elastic, 
plastic, or viscoplastic material responses, as well as thermal loading. Moreover, the domain 
integral formulation is relatively simple to implement numerically, and it is very efficient. This 
approach is very similar to the VCCT. 
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Taking account of plastic strain, thermal strain, body forces, and crack face tractions leads to the 
following general expression for / in two dimensions: 



chi; 


3xj 


3xi 


3w p dO 9uj 

+ aa x j Fj — - 

dxi J 3xj 3xj 


dA- J 

r+r 


—^q dr (2.27) 
dx 1 


where a is the coefficient of thermal expansion, <9 is the temperature relative to ambient, the 
super-script p denotes the plastic strain, and q is an arbitrary but smooth function that is equal to 
unity on r o and zero on /] (Figure 11). 

Inertia can be taken into account by incorporating T, the kinetic energy density, into the 
group of terms that are multiplied by q. For a linear or nonlinear elastic material under 
quasistatic conditions, in the absence of body forces, thermal strains, and crack face tractions, 
Equation (2.27) reduces to 



du: 

7 — 

1J 3x, 


-Mi 


dq 

dxj 


dA 


(2.28) 



Figure 1 1 . Inner and outer contours, which form a closed contour around the crack tip when connected by 

r+ and r (Anderson, 1991). 

Equation (2.28) is equivalent to Rice’s path-independent /-integral. When the sum of the 
additional terms in the more general expression (Equation (2.27)) is nonzero, / is path 
dependent. The variable, q, can be interpreted as a normalized virtual displacement. The q 
function is merely a mathematical device that enables the generation of an area integral, which is 
better suited to numerical calculations. 

Implementing /-Integral into FEA 

The domain integral approach can be implemented into finite elements with some modifications 
(Shih et al., 1986; Dodds and Vargas, 1988). In 2-D problems, one must define the area over 
which the integration is to be performed. The inner contour, r o is often taken as the crack tip, in 
which case A* corresponds to the area inside of /]. The boundary of /] should coincide with 
element boundaries. An analogous situation applies in three dimensions, where it is necessary to 
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define the volume of integration. The latter situation is somewhat more complicated, however, 
since J( r/) is usually evaluated at a number of locations along the crack front. 

In the absence of thermal strains, path-dependent plastic strains, and body forces within the 
integration volume or area, the / integral is expressed as follows: 


crack faces 


9uj 

(J 2 J— q w 


where the spatial derivatives of q are given by 


dq _y^ 2 y^ 3 dNj 
3x i £ ~ dXj 


(2.30) 


and m is the number of Gaussian points per element, w p are w are weighting factors, n is the 
number of nodes per element, q, are the nodal values of q,, ^ are the parametric coordinates for 
the element, and Ni are the element shape functions. The quantities within { } p are evaluated at 
the Gaussian points. Note that the integration over crack faces is only necessary when there are 
nonzero tractions. 

Finite Element Mesh Design 

Although many commercial codes have automatic mesh generation capabilities, producing a 
properly designed FEM always requires some human intervention. Crack problems, in particular, 
require skill on the part of the user. Figure 12 shows several common element types for crack 
problems. At the crack tip, four-sided elements (in 2-D problems) are often degenerated down to 
triangles, as Figure 13 illustrates. An analogous situation occurs in three dimensions, where a 
brick element is degenerated to a wedge. 



<b) 20 node brick element (d) 27-node brick element 


Figure 12. Isoparametric elements that are commonly used in 2-D and three-dimensional (3-D) crack 

problems. 
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Figure 13. Degeneration of a quadrilateral element into a triangle at the crack tip. 

Using the Strain Softening Method to Perform Progressive Delamination Analysis 

Using the strain-softening method, Davila et al., have developed a method to simulate 
progressive debonding or delamination using decohesion elements (Davila et al., 2001). The 
onset of damage and the growth of delamination can be simulated without previous knowledge 
about the location, size, and direction of propagation of the delaminations. The approach has 
been found to be well suited to nonlinear progressive failure analyses where both ply damage 
and delaminations are present. 

The approach consists of placing interfacial decohesion elements between composite layers. A 
decohesion failure criterion that combines aspects of strength-based analysis and fracture 
mechanics is used to simulate delamination by softening the element. As the delamination grows, 
the tractions across the interface will grow to a maximum value, then decrease and eventually 
vanish when complete decohesion occurs. The work of normal and tangential separation can be 
related to the critical values of ERR. 

Davila et al. developed an 8-node decohesion element and implemented it in ABAQUS™ to 
perform the analyses (de Moura et al., 2000; 1997). The decohesion element models the interface 
between sublaminates or between two bonded components. The element consists of a zero- 
thickness volumetric element in which the kinematics of the elements that are being connected 
on the top and bottom faces are correctly modeled. The element represents damage by a material 
response using a cohesive zone ahead of the crack tip to predict the delamination growth. The 
interface elements have been formulated to deal with mixed-mode delamination growth 
problems. 

An example comparing this approach with a beam solution and test data uses a DCB to 
determine interlaminar fracture toughness in Mode I (Gjc)- A Gr/Ep specimen made of a 
unidirectional fiber-reinforced laminate containing a thin insert at the mid-plane near the loaded 
end is chosen. The specimen was 15 cm (6 in.) long and 2 cm (0.8 in.) wide with an initial crack 
length of 5.5 cm (2.2 in.). A plot of the reaction force as a function of the applied end 
displacement, d, is shown in Figure 14. The beam solution was developed by Mi and Crisfield 
(Mi and Crisfield, 1996) for isotropic adhered materials and used plane stress assumptions. After 
the initiation of delamination, fiber bridging in the test specimen caused a small drift in the 
response compared to the FEM and analytical results. The results of the decohesion method were 
in good agreement with the test data. 
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DisoEacennenutf, mm. 

Figure 14. Load-deflection response ofDCB test. 

A parametric study revealed that the accuracy of the prediction was reduced when the element 
sizes in the region of strain softening were greater than a prescribed value. Poor results were 
obtained when the element size was greater than 1.25 mm (0.5 in.). 

2. 2. 2. 2 FEA With Singular Element Formulation 

In contrast to the conventional delamination prediction approach presented in the previous 
section, the following approach uses singular elements to predict delamination growth along a 
crack tip. 

Davidson developed a crack-tip element methodology for predicting delamination growth in 
laminated composite structures (Davidson, 1998). First, the toughness versus mode mix relation 
was determined experimentally. The definition of mode mix sought was the one that, when 
applied to all experimental results, produced a single-valued toughness versus mode mix curve. 
The mode mix definition was found from the results of a series of delamination toughness tests, 
and obtained within the construct of a crack-tip element analysis. This approach allowed for the 
possibility that the singular field-based decomposition was valid; it also allowed for the 
possibility that an alternative definition was valid. This alternative definition was based on 
parameters that uniquely defined the loading at the crack tip but which were insensitive to the 
details of the near-tip damage. The second component of the methodology consisted of analyzing 
the structure of interest and assessing delamination growth. This was also done using the crack 
tip element analysis and incorporated the experimentally determined definition of mode mix. 

One of the attractions of this approach is that it removes the need for detailed 2-D and 3-D FEAs. 

Total Energy Release Rate 

Consider the crack tip element of Figure 15 that represents a 3-D portion of the crack tip region 
in a general interfacial fracture problem. The crack lies locally in the x-y plane at constant z, the 
element is in a state of plane stress or plane strain with respect to the y coordinate direction. The 
crack front may be straight or curved with the length of the element in the y direction assumed to 
be sufficiently short such that there is no significant variation in the loading and in the 
orientation of the edges of the element with respect to y. It is assumed that the lengths of the 
cracked and uncracked regions comprising the element are large with respect to their thicknesses 
but are sufficiently small such that geometric nonlinearities are negligible. 
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Figure 15. Crack-tip element and local loading. 


Classical plate theory is used to predict the overall deformations and strain energies of the 
element (Davidson et al., 1995; Schapery and Davidson, 1990; Hu, 1995). It has been shown 
(Davidson et al., 1995; Schapery and Davidson, 1990) that the loading on the crack-tip element, 
which produces the stress singularity can be fully characterized in terms of a concentrated crack 
tip force, N c , and moment, M c . The concentrated crack-tip force and moment are found by 
enforcing the condition that the displacements of the upper and lower plates be compatible along 
the crack plane over -b < x < 0. The ERR of the crack-tip element is obtained through a 
modified virtual crack closure method and may be expressed in terms of N c and M c , rather than 
the four independent quantities Ni, N 2 , Mi, and M 2 . This gives 



( c lNc + c 2 Mg + 2^/qc7N c M c sinr ) 
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The concentrated crack-tip force and moment are given by 
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and ti and t 2 are the thicknesses of plates 1 and 2 as shown in Figure 15. 

Mode Decomposition - Classical Definition 

Mode decomposition is performed to define the ERR in terms of its Mode I and Mode II 
components. The Mode I and Mode II ERRs are given by 
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where /) is defined as 
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In the above equations, L is a fixed dimension, regardless of the body being analyzed, whereas L 
is a characteristic dimension that scales with the body. Equation (2.39) ensures that the correct 
mode mix will be predicted for those cases where an oscillatory singularity exists. That is, i2is 
determined for a specific crack-tip element geometry and is based on the characteristic 
dimension L. The scaling of Q given by Equation (2.39) ensures that the same £2 will predict 
the correct mode mix for arbitrary loadings of a geometrically similar element with different 
absolute dimensions (Davidson et al., 1995). 

Example Problems in Two Dimensions 

Instability-Related Delamination Growth 

Davidson presented the prediction of instability-related delamination growth (Davidson, 1998) 
using a global cylindrical buckling analysis along with a local crack tip element analysis 
(Davidson and Krafchak, 1993). Figure 16 shows a cross-sectional view of a midplane 
symmetric laminate containing two through-width delaminations located symmetrically with 
respect to the laminate’s midplane. 
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A cylindrical buckling analysis is performed on the laminate model (Davidson and Krafchak, 
1993; Yin, 1988). The post-buckled laminate appears as shown in Figure 17, with all internal 
force and moment resultants shown in their positive conventions (the actual in-plane forces are 
compressive in this problem and therefore the reverse of those shown in the figure). 



Delaminations 

Figure 16. Cross-sectional view of laminate containing two symmetrically located delaminations that is 

subjected to compressive loading. 



Figure 17. Post-buckled configuration. 

Figure 18 shows the crack-tip element and local loading. Minor modifications to the crack-tip 
element formulation were made to enforce the symmetry constraints. The loads, N, Ni, and N 2 , 
and moment, Mi, acting on the crack- tip element are given by 

TsjB 

N 1= N° N 2 =-^- N = Nj +N 2 M 1= M°-NxW (2.40) 

Substituting these loads into the crack-tip element equations provides an analytical methodology 
for determining mixed-mode ERRs as a function of the applied far- field loading. 

Figure 19 presents a comparison for total ERR between the linear crack-tip element analysis and 
a geometrically nonlinear finite element analysis (Davidson and Krafchak, 1993). A [02/90/02]3 S 
Gr/Ep laminate with delaminations at the interfaces between the 5 th and 6 th and 25 th and 26 th plies 
was modeled in both cases. The laminate is assumed to be in a state of plane strain. The 
horizontal axis, denoted as applied compressive strain, is the deformation of one of the loaded 
ends of the laminate divided by the laminate half-length. The value of Q used was 14.6° and 
was found by a linear finite element analysis of the crack-tip element geometry. Because the 
delamination growth occurs at a 0/0 interface, the stress singularity obeys a classical inverse- 
square root relation. 
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Figure 18. Crack-tip element and loading for delamination buckling problem. 



Figure 19. Comparison of total ERR for a [02/90/02]3 S laminate. 


Using this method, only a single, linear finite element analysis need be performed to obtain Q, 
thereafter, all results are generated analytically. For certain specific problems, no finite element 
analyses are required, as the value of Q may be taken directly from Davidson et al. (1995). 

Edge Delamination 

Consider a mid-plane symmetric laminate with a single delamination at its free edge, as shown 
in Figure 20. The term “sublaminates” will be used to refer to the regions bounded by the 
delamination and the laminate free surfaces. It is assumed that the delamination length, a, is large 
compared to the thickness of either of the sublaminates and that the strains and curvatures in the 
laminate and sublaminates are given by 

£ \o= £ o 7760 =° Zl=Zo Z6=° ( 2 - 41 ) 


First, consider the classical laminated plate theory solution to an uncracked laminate under a 
specified loading. A transverse stress distribution is obtained for which the resultant transverse 
force, N 2 , and moment, M 2 , vanish along those faces of the laminate which are defined by a 
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normal vector in the positive or negative X 2 direction. Assuming a crack is introduced at an 
arbitrary interface, the resultant transverse forces and moments along the free edge of each of the 
sublaminates, as well as the shear and normal stresses along the crack plane, must vanish. To 
achieve this, superpose onto the solution to the uncracked laminate the solution to the problem 
where the sublaminates are loaded by transverse forces and moments, which are equal and 
opposite to those which are given by the uncracked solution. This superposition problem is 
shown in Figure 21 (with all forces and moments shown in their positive conventions). 



fl» 

Figure 21. (a)Transverse stress distribution for an uncracked mid-plane symmetric laminate subjected to 
uniform axial extension in the Xi direction, (b) Crack-tip element. Superposition of problem (b) onto 

(a) gives solution to cracked laminate. 
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Since the ERR for the uncracked laminate (Figure 21(a)) is zero, the ERR for the cracked 
laminate is simply that from the “crack-tip element” of Figure 21(b). The ERR and fracture mode 
ratio for this element may be obtained through a modified virtual crack closure method 
(Davidson, 1995). The equation for G is derived as in Equation (2.31) above. The subscripted 
quantities A', B', and D' are defined such that, for sublaminate 1 or 2 

4o = AjN® + BjM®, i = 1,2 (2.42) 

= B-N® + D-M^p, i = 1,2 (2.43) 


where 

£ 2o ~ ' s midp lane strain in the x 2 direction in sublaminate 1 

X 2 * = is the curvature in the x 2 direction in sublaminate 1, etc. 

For those cases where crack growth is between orthotropic layers and the crack tip stress field 
exhibits an inverse square root singularity, G\ and Gw may each be expressed as Equations (2.37) 
and (2.38) where Q = Q. 

The value of G\ and Gw will never be less than zero. However, when the bracketed portion of G\ 
is less than zero, a negative Mode I stress intensity factor is predicted and the Mode I ERR 
component is actually a “closing mode.” This implies a crack face interference and the present 
solution, which does not enforce crack face contact constraints, is invalid for that particular 
loading case. 

The value of Q is independent of the loading and may be found from numerical or other results 
of one special loading case. This may be done, for example, by FEA of the crack-tip element or 
taken from a generalized plane strain FEA of the free edge delamination problem. Also, the 
value of the total ERR is independent of Q, that is, the sum of G\ and Gw will always give the 
same result as G, regardless of the value of Q. 

As an example, consider the determination of the ERR in a [45/0/-A5/90] s T300/5208 laminate 
with delaminations at both -45/90 interfaces (Davidson, 1994a and 1995). The laminate is 
subjected to a uniform strain in the global xi direction. The value of the ERR for various 
hygroscopic loading cases are presented in Figure 22 along with the finite element results from 
O’Brien et al. (1986). The results are found to match almost exactly. 
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Figure 22. ERR versus hygroscopic loading for a [45/0/ — 45/90] s T300/5208 laminate with edge 

delaminations at both (- 45/90). 


2.2.2.S Progressive Delamination Growth 

A self-contained finite element program developed for progressive delamination growth 
technique allows non self-similar delamination progression. The approach is demonstrated on the 
post-buckling delamination growth of a thin composite sublaminate from a thicker sublaminate 
structure shown in Figure 23. In order to save computational costs, a double shell model is 
pursued using Mindlin shell theory. Efficiency and accuracy of this approach has been 
demonstrated in comparison to a 3-D model. 



Figure 23. Configuration of embedded delamination. 

The conditions for delamination growth are as follows: (1) the delamination exists prior to 
loading, (2) the growth is along the delamination interface plane, (3) stable growth occurs by 
increasing the boundary load, and (4) the computed distribution of strain energy release rate will 
exceed the critical value at the section of the delamination front that is propagating. 


The modified crack closure technique (Rybicki and Kanninen, 1977; Jih and Sun, 1990) is used 
to calculate the strain energy release rate. In this technique, the strain energy released during 
crack extension is assumed to be equal to the work needed to close the opened surfaces. 
Generally, if the relationship between the nodal force and relative displacement is nonlinear, the 
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closure operation must be performed incrementally (Cheong and Sun, 1992). However, it was 
discovered that even while using large deflection theory, the relations between the nodal force 
(moment) and relative displacement (rotation) are nearly linear at the delamination front. Thus, 
one-step closure is adopted in this study. This condition must be verified for each problem as 
there may be problems where the linear relation does not persist. 

In the crack closure method, the virtual crack extension A a is usually taken to be small compared 
to the crack length. Thus, the nodal displacements (crack opening displacements) at the original 
crack tip after the extension can be approximated by the nodal displacements behind the crack tip 
before the extension. This reduces the analysis to one step instead of two. This procedure is valid 
if the elements in front of and directly behind the delamination crack front are similar. The 
conditions for similarity require that the size of each element both in front of and directly behind 
the crack front must be equal, the length of the elements normal to the delamination front must 
be A a, and ideally, the two sides of the elements should be made perpendicular to the 
delamination front. 

A program, FECM209, developed at Purdue University is used to simulate delamination 
propagation. For simplicity, a total strain-energy release rate, G, is chosen for the approach. The 
procedure for progressive delamination growth is as follows: 

1 . Calculate the linear stiffness matrix. 

2. Apply the in-plane load incrementally by a specified amount Aeo until the final load is 
reached. 

3. Compute the nonlinear stiffness matrix based on the current stress field and 
displacements. 

4. Solve the displacement and stress fields. 

5. Find the converged solution or return to Step 3 if the solution is not convergent. 

6. Evaluate the strain-energy release rate along the delamination front. 

7. Check the distribution of the strain-energy release rate to find the comer nodes that 
exceed the critical value. If none exists, return to Step 2. 

8. Move the nodes in the outward normal direction to the delamination front where the 
critical strain-energy release rate is exceeded and generate a new mesh. 

9. Compute the linear stiffness matrix for the current mesh. 

10. Go to Step 3. 

In Step 6 the local value of G is computed at every comer node. When the strain-energy release 
rate at a node exceeds G cr , the node will propagate in the local outward normal directions. The 
amount of propagation is proportional to the difference in the strain energy release rate and G cr . 

During remeshing, the ratios of element sizes along the delamination front remain constant. In 
addition, the number of elements along the x and y axes remains constant. A large change in 
element sizing could lead to poor results. 

One problem with changing the mesh during loading is that the stiffness matrix must be updated 
accordingly. In Step 9, the linear stiffness matrix is updated based on the current mesh, but the 
nonlinear stiffness matrix must be calculated based on current displacements and stress fields. 
Since the location of the nodes has changed, the displacements and stresses at the new location 
should be interpolated from those at the previous locations. However, since the movements of 
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the nodes are very small each time the mesh is updated, the displacements and stresses at the 
previous location can be used as an approximation for the current location. 

Most commonly used commercial codes do not allow the movement of nodes in a mesh during 
the loading procedure. Consequently, after the mesh is updated, the loading procedure would 
have to be restarted. In other words, the whole loading process has to be restarted for each 
movement of the delamination front. This is computationally very time-consuming for 
simulating delamination growth. 

Post-buckling Delamination Growth Predictions 


The following are results from simulations of post-buckling embedded and edge delaminations 
using the above technique. An example of an embedded circular delamination is shown in 
Figure 24. An initially circular delamination front with a radius of 15 mm (0.6 in.) is 
embedded in a 100 mm (3.9 in.) square panel. The critical strain-energy release rate 
G cr = 200 J/m 2 (13.7 ft-lb/ft 2 )is chosen for the analysis. The panel is loaded in the x-direction 
under uniform displacement. 

Overlapping of the upper and lower sublaminates often occurs over some range at the 
delamination front. This overlapping can be constrained by contact elements. 


The growth contours and corresponding normalized strain-energy release rate distributions for 
the [(90/0) s ] sublaminate are shown in Figures 24 and 25, respectively. The upper sublaminate is 
four plies while the lower sublaminate is twenty-eight plies. Transverse deflection of the lower 
sublaminate is constrained. The growth initiates perpendicular to the loading direction and is 
stable for the first four loading levels (£\ - £4). Curves labeled £5 and £4 have become neutrally 
stable since the loading level is not increased and the distribution of the strain-energy release rate 
in the growth-region remains near G cx . The last curve labeled £1 is unstable because G has 
increased to 208 J/m 2 (13.8 ft-lb/ft 2 ) while the loading level remained constant. 
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Figure 24. Progressive delamination fronts for the 
[(90/0) s ] laminate. 



Figure 25. Strain-energy release rate distributions 
for the [(90/0) s ] laminate. 


The distribution plot of G for the ideal ellipse (with the aspect ratio of contour 86) has a very 
large oscillation about the G cr line while the distribution at £e is constant. In fact, the maximum 
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value is more than twice the value of G a . Thus, the propagating delamination cannot be modeled 
as an elliptical crack. 

It has been found that the initiation and subsequent stable delamination growth is very sensitive 
to the stacking sequence of the upper sublaminate. Lower bending stiffness parallel to the 
loading direction leads to higher predictions of strain-energy release rate and delamination 
growth at lower load levels. Un-symmetric sublaminates may have strong extensional-bending 
coupling that effectively reduces the bending stiffness leading to higher strain energy release 
rates. 

The analysis neglects matrix cracking that can occur in the post-buckled sublaminate as observed 
by Whitcomb (1988). Matrix cracking can change the stiffness and may affect the stability of the 
growth. 

A similar approach has been used to simulate the damage at the tip of a notch in cross-ply 
laminates. A schematic of the failure mode is shown in Figures 26 and 27. An axial split forms in 
the 0°-ply from a matrix crack and a delamination forms between the 0° and 90° layers. The 
fractures at the split and along the delamination front are modeled independently. This model has 
been used to investigate the size of the damage zone under various degrees of ply clustering or 
repeated plies. Larger damage zones are predicted and observed for increased ply clustering. 
Laminate strength approaches the laminate net section strength in the with large-ply clustering 
and damage zones. 

See “Efficient Modeling of Postbuckling Delamination Growth in Composite Laminates Using 
Plate Elements” by Klug et al. (1996) for more detailed descriptions of the analysis results. 



Figure 26. Schematic of axial splitting and delamination in a [(90/0)] s laminate. 
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Figure 27. Local finite element mesh. 


2.3 Impact Damage 

2.3.1 Prediction of Burst Strength After Impact (BAI) for Composite Overwrapped Pressure 
Vessel (COPV) 

The prediction of BAI for COPVs is one of most challenging tasks in COPV fracture control. 
Chen et al. (1998 and 1999) developed a procedure that can be used in the BAI predictions for 
COPVs under impact. A 3-D FEM was constructed with multi-layered, laminated, composite 
shell elements for the cylindrical COPV and the analysis was performed with ABAQUS™ 
(Vaidya et al, 1998). The model consisted of two portions. In the cylindrical portion, where the 
low energy impact occurs, the liner thickness and the composite overwrap lay-up were used to 
obtain individual ply stress/strain responses. In the hemispherical ends, where the laminate 
material thickness and fiber orientation change continuously, integrated properties for the 
overwrap were used for simplicity. Due to symmetry, only half of the vessel was modeled, as 
shown in Figure 28. The vessel was fully constrained on one end while the other end was 
constrained only in the radial and tangential directions to simulate the constraint condition in the 
tests. 





Figure 28. FEM of cylindrical COPV under impact. 
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The impactor mass was specified through the use of a point mass element, and the impactor 
geometry was specified using a rigid surface element. Surface definitions were specified for both 
the impactor and the COPV to ensure proper surface interaction during impact. The analysis was 
carried out in two steps. In the first step, an internal pressure was applied to the shell structure, 
and its static response was determined. In the second step, the impactor was assigned an initial 
velocity, and the resulting transient dynamic response of the COPV to the impact was evaluated 
using the implicit time integration procedure provided in ABAQUS™. 

2. 3. 1.1 Impact Damage Prediction 

Since fiber breakage is more detrimental to the tensile fracture strength of composite laminates 
than either delamination or transverse matrix cracking, different modes of impact damage must 
be distinguished in order to accurately predict the BAI of an impact damaged COPV. A modified 
maximum strain failure theory was used to predict both the extent and the mode of composite 
damage. The ply normal strains, in the fiber and transverse directions, and shear strains at each 
integration point were calculated and compared with their respective ultimate strains to 
determine whether any damage occurred in each ply. The three distinctive damage modes were: 
damage in the fiber direction (Mode I), damage in the transverse direction (Mode II), and 
damage in shear (Mode S). The extent of damage was determined by examining the strains in all 
the finite elements adjacent to the impact location. 

A FEM of the cylindrical COPV is shown in Figure 28. The predicted damage mode and extent 
of damage obtained from the FEA results are summarized in Table 1, where the inner ply is 
shown at the top of the table and outer ply is shown at the bottom. The damage modes for the 
two elements adjacent to the impact location are separated by a dash, and a “0” indicates that 
there is no damage. A blank space in the table indicates that neither element was damaged and a 
“2s” means that both Mode II and Mode S damages were predicted. It is observed that the 
damage concentrates on the 72.5° plies and the 27.7° plies. For the 15 and 20 ft-lb (20 and 
27 joules) impacts, the inner 19.2° ply was also damaged. 


Table 1. 

Predicted Damage Modes for Each Ply 


Impact Energy 

Laminate 

5 ft-lb 

10 ft-lb 

15 ft-lb 

20 ft-lb 

Ply 

(7 joule) 

(14 joule) 

(20 joule) 

(27 joule) 

19.2° 



0-S 

0-S 

-19.2° 





85.7° 





-85.7° 





27.7° 



1-0 

1-S 

-27.7° 


S-0 

S-1 

S-1 

72.5° 


S-2 

2S-2 

2S-2 

-72.5° 

2-0 

2-2S 

2-2S 

2-2S 


85.8° 

-85.8° 


To account for the effect of damage on laminate strength, empirical material degradation factors 
(DF) were used to reduce the stiffness of damaged plies in the subsequent strength calculations. 
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The selection of the DF was based on the following empirical rules. For undamaged plies, the DF 
was zero. For plies with damage, the DF is as follows: 

DF = 1.0 if damaged in Mode I 
DF = 0.8 if damaged in Mode I or IIS 
DF = 0.6 if damaged in Mode S 

Mode IIS is considered to be no worse than Mode II and therefore the same DF was assigned to 
these two modes. In this way, the material effectiveness can be obtained for each ply by 
considering the damage mode and the damage extent, as 

EF = 1.0-]T^ (2.44) 

where EF is the ply material effectiveness factor and summation is performed over the two 
elements surrounding the impact location. 

For the COPV studied, the results of the material effectiveness calculation of each ply are 
presented in Table 2. It is observed that the material effectiveness in one of the 72.5° plies was 
reduced to 60 percent under the 5 ft-lb (7 joule) impact. Similarly, the 72.5° plies and the 27.7° 
plies were all reduced to 20 percent effective, and the inner 19.2° ply was reduced to 70 percent 
effective under the 20 ft-lb (27 joule) impact. 

Table 2. Calculated Ply Material Effectiveness Factors 


Laminate 

Ply 

Impact Energy 

5 ft-lb 
(7 joule) 

10 ft-lb 
(14 joule) 

15 ft-lb 
(20 joule) 

20 ft-lb 
(27 joule) 

19 . 2 ° 

- 

- 

0.7 

0.7 

- 19 . 2 ° 

- 

- 

- 

- 

85 . 7 ° 

- 

- 

- 

- 

- 85 . 7 ° 

- 

- 

- 

- 

27 . 7 ° 

- 

- 

0.5 

0.2 

- 27 . 7 ° 

- 

0.7 

0.2 

0.2 

72 . 5 ° 

- 

0.3 

0.2 

0.2 

- 72 . 5 ° 

0.6 

0.2 

0.2 

0.2 

85 . 8 ° 

- 

- 

- 

- 

- 85 . 8 ° 

- 

- 

- 

- 


2. 3. 1.2 Post Impact Burst Strength Predictions 

It was observed in the impact-testing program that debonding between the aluminum liner and 
the composite laminate occurred upon impact. Therefore, the aluminum liner and composite 
laminate were assumed to be un-bonded in the analysis and their contributions to the BAI of the 
COPV were calculated separately. The BAI of the COPV is the sum of the two component 
strengths. Once the impact damage was determined, the BAI of the COPV laminate was 
determined by a progressive failure analysis procedure. As the aluminum liner has yielded long 
before the laminate fails, the burst strength contribution by the aluminum liner, p, was 
conservatively estimated with a closed-form formula 



r 
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(2.45) 



where r is the radius of the pressure vessel; s y and t are the aluminum liner’s yield strength and 
thickness, respectively. The yield strength of the aluminum liner was assumed to be 42,000 psi 
(290 MPa). 

The burst strength contributed by the undamaged laminate was predicted to be 10,170 psi 
(70 MPa). The contribution of the aluminum liner was approximately 510 psi (3.5 MPa). The 
total predicted burst strength was therefore 10,680 psi (73.6 MPa), which was in close agreement 
with the average test value of 10,700 psi (73.8 MPa). Table 3 compares the predicted burst 
strengths with test data for the cylindrical COPV with damage caused by impact at various 
energy levels. The degradation in the burst strength was predicted to range from 7 percent to 
29 percent for impact energies of 5 to 20 ft-lb (7 to 27 joules), respectively. One of the 15 ft-lb 
(20 joule) BAls does not appear to be reasonable because it indicates a greater residual burst 
strength than that of the 10 ft-lb test (14 joule). Ignoring that particular test data, one can see that 
the predicted degradation in burst strength was a little larger than the test data for the 10 ft-lb 
(14 joule) case and fairly close to the test data for all other cases. (See “Impact Damage Effects 
on Gr/Ep Composite Overwrapped Pressure Vessels” by Chen et al. (1999) for more detailed 
descriptions of the test results.) 

Table 3. Comparison of Burst Strengths for Cylindrical COPV 


Impact 

Test Results 

Analytical Results 

Energy 

Burst Pressure 

Degradation 

Burst Pressure 

Degradation 

0 ft-lb 

10,700 psi 


10,680 psi 


(0 joule) 

(73.7 MPa) 


(73.6 MPa) 


5 ft-lb 

9,800 psi 

8.4% 

9,920 psi 

7.1% 

(7 joule) 

(67.6 MPa) 

(68.4 MPa) 

10 ft-lb 
(14 joule) 

8,884 psi 

17.0% 

8,220 psi 

23.0% 

(61.2 MPa) 

(56.7 MPa) 

15 ft-lb 

8,246 psi 

22.9% 

7,960 psi 

25.5% 

(20 joule) 

(56.9 MPa) 

(54.9 MPa) 

15 ft-lb 

8,377 psi 

21.7% 

7,960 psi 

25.5% 

(20 joule) 

(57.8 MPa) 

(54.9 MPa) 

15 ft-lb 

9,257 psi 

13.5% 

7,960 psi 

25.5% 

(20 joule) 

(63.8 MPa) 

(54.9 MPa) 

20 ft-lb 

7,681 psi 

28.2% 

7,590 psi 

28.9% 

(27 joule) 

(53.0 MPa) 

(52.3 MPa) 


2.3.2 Predicting Damage after Low- Velocity Impact 

Choi and Chang developed a method for predicting the impact damage of laminated composite 
plates resulting from both line-nose and point-nose impacts (ABAQUS™ User’s Manual, 1992). 
The method consists of a stress analysis for determining the stress distributions inside the 
laminates during impact and a failure analysis for predicting the initiation and extent of the 
impact damage. The authors used a transient dynamic stress analysis code developed by Wu and 
Chang (Choi and Chang, 1992). The stress analysis was based on a 3-D linear elastic theory that 
assumes the materials in each layer to be homogeneous and orthotropic. Eight-node brick 
elements incorporating incompatible modes were used in the FEA. Elements with incompatible 
modes were chosen because of the improved accuracy in calculating bending stiffnesses and 
interlaminar shear stresses. The method was found to correspond well with experiments in 
predicting the location and extent of impact damage. 
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The procedure for determining the extent of the impact damage was as follows: 

1 . Calculate the transient dynamic stresses within each layer as a function of time. 

2. Apply the matrix failure criterion for predicting the critical matrix cracks in each layer for 
determination of the extent of delaminations. 

3. If matrix cracking is predicted in a layer of the laminate, then apply the delamination 
criterion subsequently in the upper and lower layer of the interface during the entire 
period of impact. 

The procedure should be repeated at the other layers during impact for determining any 
additional matrix cracking and delaminations. The final size of each delamination is determined 
by the area within which the stress components satisfy the delamination failure criterion during 
the entire duration of the impact. No material degradation was considered in the model, and the 
model does not take into account delamination interaction during impact. 

Two failure criteria were proposed in the study to predict the initiation of critical matrix cracks 
and estimate the extent of delamination. 


2.3.2. 1 Critical Matrix Cracking Criterion 

The failure criterion proposed for matrix cracking can be expressed as (Wu and Chang, 1989) 

fe M >l Failure 


n = ) 

1 

fn= 3 

<7yy 

+ 

^yz 

n Y 


n Si 

J 


v 1 J 


:e M 


e M < 1 No Failure 
n Y= n Y t if d\,„ >0 


(2.46) 
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n Y= n Y c ifdyy <0 


where x and y are local coordinates of the nth layer parallel and normal to the fiber directions, 
respectively, and z is the out-of-plane direction. S, is the in situ interlaminar shear strength within 
the laminate under consideration and Y t and Y c are the in situ ply transverse tensile and 
compressive strengths, respectively, as defined by Choi and Chang (ABAQUS™ user’s manual, 
1992). 

Whenever the calculated average stresses in any one of the plies in the laminate first satisfy the 
criterion (i.e., eM = 1) during impact, the initial impact damage is predicted. It was assumed that 
the matrix crack would propagate throughout the thickness of the ply group that contained the 
cracked ply. A delamination can be immediately induced from the location of the matrix crack 
along the interfaces of the ply group. Because additional matrix cracking can occur in other 
layers as time increases, the criterion should be continuously applied to the neighboring layers. 

2. 3. 2. 2 Impact-Induced Delamination Criterion 

Once the crack matrix criterion is applied and a critical matrix crack is predicted within a 
laminate, a delamination is initiated from the crack. A semi-empirical model was developed to 
estimate the extent of delaminations in the composites. The model includes the following two 
types of critical cracks initiating delaminations resulting from the impact: the shear crack 
generated within the laminates, and the bending crack produced from the bottom surface of the 
laminates as shown in Figures 29 and 30. 
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It was reasoned that delamination growth due to low- velocity impact could occur only when two 
sequential conditions were met by the following: 

1 . One of the ply groups directly above or below the concerned interface failed due to 
matrix cracking. 

2. The combined stresses governing the delamination growth mechanisms through the 
thicknesses of the upper and lower ply groups of the interface reaches a critical value. 
This is the impact-induced delamination growth criterion. 



Figure 29. Basic impact damage mechanism of laminates subjected to 2-D line-loading impact. Top: 
delamination induced by inner shear cracks. Bottom: delamination induced by surface bending crack. 


Impact Damage Growth Mechanism 



Figure 30. Schematics of the delamination growth mechanism due to a shear crack and a bending crack, 
respectively, in a laminated composite subjected to point-nose impact. 


The second condition can be written in equation form as 
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(2.47) 


where D a is an empirical constant which must be determined experimentally. The constant has 
been found to be insensitive to ply orientation and thickness of the laminate, and primarily 
dependent on the material system used. The superscripts n and n + 1 correspond to the upper and 
lower plies of the nth interface, respectively. er yz and <Xyy are the averaged interlaminar and in- 
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plane transverse stresses within the nth and n + 1th ply, respectively. <j xz is the averaged 
interlaminar longitudinal stress within the n + 1th ply. 

3.0 Damage Tolerance Life Prediction Methodologies 

This section presents life prediction methods for damage growth under fatigue-type loading. 

3.1 Delamination Growth in Cyclic Loading 

3.1.1 Closed Form Solution for Delamination Growth 

A delamination growth criteria that incorporates Mode II was investigated (dL/dN = f(Gn)). A 
different delamination growth-rate criterion has been proposed that includes the effects of both 
A Gii and Gii max . The following expression was found to provide best fit of the experimental data 
for a wide range of stress ratios: 



(3.1) 


This equation is similar to the approaches used to characterize crack growth rates in metals 
for different stress ratios (Broek, 1986), and for R = °° it reduces to Equation (3.1) since 
Gumax = AGii. A comparison of the analytical and experimental cyclic growth rates for three 
different R values is shown in Figure 31. Again, good agreement in the cyclic growth rates is 
observed for larger delaminations for which Mode II drives damage growth. 



,L [mm] 

Figure 31. Cyclic growth rates (dL/dN) in constant-amplitude compression-compression loading as a 
function of L (solid lines - experimental data, dashed lines - analytical prediction). 

A 3-D FEM would be required for a full understanding of parameters that influence delamination 
growth. However, the simplified modeling approach given by Han et al. (1999) can provide 
useful information about the major forces that drive delamination growth. For a delamination 
that spans only a few surface plies the results suggest that Mode II is the dominant factor that 
governs delamination growth with minor contributions from Mode I for smaller delamination 
lengths. If the delamination is nested deep in the thickness of the laminate, the influence of 
Mode I component becomes dominant and cannot be neglected (see Figure 32). However, the 
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maximum value of G\ (Whitcomb, 1981). Furthermore, large scatter in reported critical strain 
ERRs (Shyprykevich, 1997) and uncertainties in loading spectrum and damage characteristics 
could lead to substantially different delamination growth predictions. 
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Figure 32. Strain ERRs (G\ and <7 n ) as a function of applied load ( (\ in percent) and delamination length 

(L) for delamination located between plies 3 and 4. 

3.1.2 Delamination Growth Using FEA 

While a number of researchers have investigated delamination growth using the concepts of 
LEFM (Whitcomb, 1989; 1990; 1992; O’Brien et al., 1986; Davidson, 1994b; Kardomateas 
et al., 1994; Russell and Street, 1989; and Tratt, 1991), the question of characterizing 
delamination growth under fatigue is still not well understood. Typically, delamination growth is 
governed by strain ERRs GI, GII, and GUI, which in the Han et al. study were determined by 
FEA (Han et al., 1999). 

Initial impact damage usually has a circular shape with multiple delaminations throughout the 
thickness. However Han argued that damage growth in fatigue more closely resembles a single 
delamination propagation front, rather than the growth of embedded multiple delaminations, as 
shown in Figure 33. Thus, an idealized plane-strain FEM was used to determine Mode I and 
Mode II components of the strain ERR. 
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Figure 33. FEM procedure. 


When subjected to a compressive load, the thin sublaminate buckles causing high interlaminar 
stresses at the delamination tip causing delamination growth. Hence, the influence of buckling 
and postbuckling of the sublaminate on delamination growth can be studied using a 
geometrically nonlinear static analysis. To initiate out-of-plane displacements an initial 
imperfection was assumed 


2a 


(3.2) 


where the maximum initial imperfection (w max ) is chosen to be 10 percent of the delaminated ply 
thickness. Strain ERRs along the delamination front were determined using a modified crack 
closure technique (Rybicki and Kanninen 1977), and the individual strain ERRs, associated with 
the two modes of fracture, were (Mukherjee et al., 1994). 
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3. 1.2.1 Comparison with Experimental Results 

Experimental results indicate that once the damage propagated, it resembled a single 
delamination growth located one or two plies below the laminate surface. To gain a physical 
understanding of the mechanisms contributing to the delamination propagation, the correlation 
between the cyclic growth rates and the strain ERRs (G) was investigated. The G\ and G\\ values 
of Figures 34 and 35 were calculated at different strain levels for delaminations located at one 
ply and two plies beneath the surface, respectively. 

In Figure 34, for a through width delamination, G\ values of a single ply are very small due to the 
low bending stiffness of the delaminated ply. As the thickness of the delaminated ply increases, 
G\ also increases (Whitcomb, 1981)as shown in Figure 35. However, the influence of Mode I 
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decreases rapidly as crack length increases, and it does not contribute to the delamination growth 
once the delamination reaches ~18 mm (0.7 inches). On the other hand, values of Gw remain 
relatively constant over the range of delamination lengths investigated; therefore, a constant 
value of Gu can be assumed for a specific value of applied strain. However, after the initial 
buckling of the delaminated plies, the crack tip actually closes due to the compressive forces that 
develop near the crack tip. These compressive forces might decrease the calculated values for 
Ci i, thus leading to the decay of cyclic growth rate. 



Figure 34. Strain ERRs (G H ) as a function of applied load (£ x in %) and delamination length (L) for 

delamination located between plies 1 and 2. 



a) Mode I. b) Mode II. 

Figure 35. Strain ERRs as a function of applied load ( t\ in percent) and delamination length (L) for 

delamination between plies 2 and 3. 


3.2 Damage Growth at Notch Tip 

In a study to predict damage growth at a notch tip of a center-notched carbon fiber/epoxy 
laminates, Spearing et al. modeled the damage as a series of interacting matrix cracks in various 
forms: splitting, delamination and transverse play cracking (Spearing et al., 1992). The basis of 
the growth model is the fatigue crack growth law known as Paris’ law (Hertzberg and Manson, 
1980). 


— = MAtf) m (3.4) 

dN 1V ’ 

where da/dN is the crack growth rate, A K is the stress intensity factor range and A.| and m are 
empirical constants. 
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Fatigue damage at notches can be characterized by the split length, 1, at the notch tip. It has been 
observed that splits and delaminations grow in combination at the notch tip. As the split length 
increases, the associated delamination area grows with a dependence on the square of the split 
length (Figure 36) which implies an increasing resistance to further crack advance. Because of 
this, Equation (3.5) can be written as 


A 

dN 


-^3 


AG 

G c 


l m /2 


(3.5) 


where 1 is the split length at the notch tip, AC is the driving force, C c the current toughness (as 
measured in monotonic loading), and A .3 is a constant. 



Split Surface 



Figure 36. Schematic of the notch tip damage pattern in a cross-ply specimen. 


3.2.1 Damage Growth in Monotonic Loading 

The model that Spearing et al. (1992) apply to damage is based on a damage growth model for 
matrix cracks in monotonic loading. Ignoring transverse ply cracking and the effect of residual 
stresses, damage at the notch in a [90/0] s laminate will grow by the simultaneous formation of 
splits in the 0° plies and triangular-shaped delamination zones at the [90/0] interfaces. 

Spearing et al. used a modified FEM developed by Kortschot and Beaumont in order to calculate 
dC/dl, the ratio of specimen compliance to an increment of split growth (Kortschot and 
Beaumont, 1991). The Kortschot and Beaumont model consists of two superimposed identical 
layers of 2-D plane stress elements. The two layers share the same nodes everywhere except in 
the delaminated region, where separate nodes are generated in the 90° ply. This model was 
further discussed in section 2. 1.2.1. The model was first developed for edge-notched specimens. 
In order to use this model for center-notched specimens, some modifications were made. The 
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relationship between the value of 3C/31 obtained from the FEM (3C/31 |fe) and those for the 
actual [90/0] s specimen is given by the scaling equation 


3C_ (W/2) FE (2t) FE 3c 
31 (W/2)(2t) 31 


where W is width and t is the thickness of the split. 

This can be extended to all laminates of the form [90i/0i] ns in which there is an equal total 
thickness of 90° and 0° plies. The scaling equation, in this case, becomes 

3C _ (W/2) fe (int) FE })C 
31 “ (W/2)(2int) 31 FE 


The method can be extended to laminates with unequal thicknesses of 0° and 90° plies by altering 
the relative ply thicknesses of the layers of the finite element mesh. 

3.2.2 Damage Growth Model Applied to [90/0] s Laminates 

For an initial length, lo, the split length after N cycles of constant load amplitude is given by 
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where G s and Gd are the energies absorbed per unit area of split and delamination, respectively, 
and AC is given by 


A G = 


(AP) 2 3C 
2t 31 


(3.9) 


where P is the range of applied load. 

Figure 37 shows a comparison of experimental data with the results from Equation (3.9) for a 
[90/0] s laminate cycled at three different peak stresses. Good agreement is shown with the 
experimental data. 



Figure 37. Normalized split length for a [90/0]s laminate as a function of load cycles showing the effect 
of varying peak stress. The solid lines represent the prediction of Equation (6.9). 
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3.2.3 Damage Growth in [90i/0i] s Laminates 


For laminates containing plies of equal thicknesses, AG increases linearly with ply thickness; 
i x t. For the [90/0] s case, there are only two interfaces at which delaminations can grow. In this 
case, the ratio AG!G C increases as i increases; so the fatigue law becomes 
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(3.10) 


Good agreement has also been found with experimental data for this case as is shown in 
Figure 38. 



Figure 38. Normalized split length as a function of load cycles for [90j/0j] ns laminates. The solid lines 

represent Equation (3.11). 


3.2.4 Fatigue Damage Growth in [90i/0j] s Laminates 


Laminates with plies of unequal thicknesses can be modeled in a similar way. The relative ply 
thicknesses of the layers in the finite element meshes were adjusted and the resulting values of 
3C/31|fe were scaled appropriately. In this case, the fatigue law is 
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From this law, it is shown that for laminates with j > i more rapid damage progression occurs 
than in [90/0] s laminates. 

3.2.5 Damage Growth in [90/0] ns Laminates 

In the case of [90/0] ns laminates, ply thickness is constant and equal for the 90° and OGn laminae; 
only the number of 90/0 interfaces varies. Through half the laminate thickness, 2n - 1 interfaces 
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exist at which delaminations can propagate. Since the split area is n x t x 1, and 3C/91 scales with 
the number of plies as shown in Equation (3.8), the fatigue law becomes 


dN 


t(Ao.) 2 (W/2)(2nt)(W/2) pE (2t) pE ^ 


-im/2 


FE 
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(3.12) 


As the number of plies increases, the model predicts that the rate of damage propagation 
decreases slightly which causes the model to underestimate the damage growth (Figure 39). This 
effect was theorized to be due to unequal damage growth through the thickness of the laminate. 



Figure 39. Normalized split length as a function of load cycles for [90/0] ns laminates. The solid lines 

represent Equation (3.13). 


3.2.6 Damage Growth in [90/+45/-45/0] s Laminates 


For these laminates, it was assumed that damage is controlled by splitting in the 0° plies and by 
delamination at the interface between the 0° ply and the innermost -45° plies. In this case, the 
finite element representation provides values for dC/dl that are independent of the split length, 
which is different from the cross-ply laminates discussed previously. The fatigue law was 
derived for this configuration in a similar manner to that of the [90/0] s laminate, yielding 
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where k is the dimension defining the intersection of the delamination with the center line of the 
notch. 

A good correlation between the experimental data and the fatigue law were also found in this 
case, as shown in Figure 40. 
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Figure 40. Normalized split length as a function of load cycles in [90/+45/-45/0] s laminates. The solid 

lines represent Equation (3.14) in an integrated form. 


4.0 Damage Tolerance Analysis Codes 

Two computer codes have been developed at The Aerospace Corporation (Aerospace) for 
predicting residual strength of composite structures containing defects. The first such code is 
Aerospace Damage Analysis Methodology for Composites (ADAM-C). The other code is 
Progressive Failure Analysis (PFA). The following sections briefly introduce these two codes. 
Detailed discussions are presented in Appendices A and B. 

4.1 ADAM-C 

In a research and development project conducted by Aerospace in the 1980s, a few analytical 
models were selected for inclusion in a static residual strength prediction software code, 
ADAM-C (Chang et al., 1985); see Figure 41 for program flow. Refer to Appendix A for a more 
detailed explanation of the study. Table 4 contains a summary of some of the analytical models 
reviewed. 



Figure 41. ADAM-C program flow chart. 
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Table 4. Assessment of Analytical Methods 


Advantages 


Disadvantages 


Loading and Damage 
Type 







Analytical Method 

Point stress failure model 

Notched fracture stress 
model 

Damage zone model 


Inherent flaw model 

Strain energy density fac- 
tor method 

Finite element approach 
Damage zone model 


Finite element method 

Edge delamination model 
Finite element method 

Thin-film 1-D model 
Finite element method 

Thin-film 2-D delamination 
model 


Simple calculation, ade- 
quate database 


Provides more accurate 
predictions 


Simple calculation, ade- 
quate database 


Provides more accurate 
predictions 

Provides adequate accu- 
racy for complex 
geometries 

Provides more accurate 
predictions 


Provides adequate accu- 
racy for complex 
geometries 


Provides adequate accu- 
racy for complex 
geometries 


Provides adequate accu- 
racy for complex 
geometries 


Stress intensity zone size is 
independent of hole sizes 

Extensive testing, limited 
database 

Lengthy calculation, not enough 
database 


Characteristic zone size is inde- 
pendent of crack sizes 


Very limited database 


Length calculation 


Prediction accuracy needs to be 
improved 


Simple calculation 


Simple calculation 


Length calculation 


Length calculation 


Simple calculation 


Length calculation 


Limited application 


Length calculation 
Limited application 


Simple calculation 


The IFM developed by Waddoups et al. was selected for use in predicting the residual strength 
of cracked composite laminates (Waddoups et al., 1971). The PSM developed by Whitney and 
Nuismer was selected for use with laminates containing holes and cutouts (Whitney and 
Nuismer, 1974). These models require simple calculations and have the advantage of an 
extensive test database for the more common composite systems such as graphite/epoxy (Gr/Ep). 
Details of these methods can be found in section 5. 1.1.1 of Whitney and Nuismer, 1974. 

For composites with delaminations at free edges, the edge delamination model (O’Brien, 1982) 
of O’Brien was selected. Through-the-width or embedded de laminations are characterized by 
Chai’s thin fdm model (Chai et al. 1981). These two models have two common features in that 
they avoid lengthy finite element calculations, and use the strain ERR G to characterize 
delamination growth. 
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4.2 PFA for Composites 

In recent years, the research and development for the PFA of composites at Aerospace has 
integrated and extended the PFA developed at NASA Langley Research Center (LaRC) in the 
late 1990s and early 2000s. The PFA is implemented in a software tool that integrates with 
existing commercial FEM software. The tool simulates the growth and propagation of failure 
modes in a composite and it is constantly adapted for new classes of problems and extended to 
incorporate new physics. The key features of the PFA code include: (1) delamination simulation 
by using interface elements; (2) Prediction of in-plane and out-of-plane nonlinear shear damage; 
(3) Prediction of matrix-cracks based on a modified version of Puck’s criterion (Puck and 
Schurmann, 1998); and (4) A fiber failure criterion based on the work of Chang and Chang 
(Chang, and Chang, 1987). 

The PFA code utilizes concepts from damage mechanics to degrade the material properties of the 
composite appropriately. The code also handles contact in the event of delamination closure 
following growth. Similar to ADAM-C, the code was developed to be as generic as possible and 
applicable to a large range of problems with adequate accuracy and conservatism. The accuracy 
of the code, however, relies on developing a complete mechanical, strength, and fracture 
database that can be used directly as input to the code. Appendix B discusses the progressive 
failure philosophy, the damage models in the PFA code, the characterization testing that is 
required to define the damage model, and validation of the PFA code. The PFA code is also 
demonstrated for few examples pertaining to pressure vessels. 
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Appendix A: ADAM-C Computer Code 
A-l ADAM-C Methodology 

Aerospace developed a computer software, ADAM-C, in the mid-1980s (Ref. A.l) with the 
intention to predict residual strength for composite structures containing defects such as cuts, 
holes and delaminations. A detailed discussions of the models incorporated in ADAM-C are 
presented as follows: 

Four basic failure modes determine the damage tolerance of composite materials: matrix 
cracking, fiber breakage, interfacial debonding, and interply delamination. All failure modes may 
be traced to localized flaws and defects that are introduced into the material through many 
potential sources, e.g., manufacturing and assembly procedures, impact damage, and 
hygrothermal environment. 

The literature survey results indicate that the growth behavior and failure modes of all flaws are 
strongly dependent on such parameters as laminate strength, stiffness, ply orientation, stacking 
sequence, fiber volume fractions, load application, and environmental conditions. These 
parameters must be integrated into a consistent methodology before full exploitation of 
composite materials can be achieved with confidence. 

A survey of aircraft industry practices revealed that the damage tolerance capability of critical 
composite hardware is often demonstrated by conducting safe-life tests on coupons, or small 
structural components, with preconceived and prefabricated flaws. This approach, however, may 
require an extensive test program when several potential laminates with different layups are 
being considered. Inasmuch as each laminate has unique properties, the test results for one layup 
are not directly applicable to the other candidate laminates. 

Several investigators have developed analytical models that predict laminate response to a 
variety of damage types common to aircraft structures. The environmental and structural 
requirements of aircraft and spacecraft, however, are quite different. Spacecraft components 
generally tend to be more highly stressed than their aircraft counterparts; hence, the applicability 
of these models to space hardware needs further evaluation. 

It is common practice in the spacecraft industry to perform acceptance proof tests on all primary 
composite structures, an approach that is both costly and time consuming. In some cases, such as 
the cryogenic temperature environment, only limited testing is possible because of the difficulty 
in simulating environmental conditions. Thus, the need for a reliable damage tolerance 
methodology that could minimize the number of these expensive tests is apparent. 

As shown in Table A-l survey results have yielded two broad categories that encompass most 
combinations of flaw type and loading condition. The category I models concentrate on through- 
the -thickness flaws such as holes, cutouts, and cracks, whereas the category II models are 
restricted to delamination, which is the predominant failure mode unique to composite structures. 
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Table A-l. Summary of Literature Survey Results 


Category 

Method/Model 

Description 

Developer 

Ref. 

No. 

1 

Inherent flaw model 

Uses linear fracture mechanics parameter for predicting 
residual strength of cracked body. 

Waddoup et 
al., 1971 

A.2 

Point stress failure 
model 

Uses laminate theory to calculate stress concentration 
factor for predicting residual strength of holes and 
cutouts. 

Whitney and 
Nuismer, 1974 

A. 3 

Finite element 
method 

Uses finite element analysis for predicting residual 
strength of laminate containing cracks. 

Wang, 1974 

A. 4 

Notched fracture 
stress model 

Relates fracture stress to composite fracture toughness; 
for cracked body residual strength calculation. 

Lin, 1977 

A. 5 

Strain energy den- 
sity factor method 

Uses strain energy density factor to predict residual 
strength of cracked composite structures. 

Sih, 1978 

A. 6 

Strain failure 
criterion 

Relates critical fiber strain to a general fracture toughness 
parameter; for cracked body residual strength prediction. 

Poe and 
Sova, 1980 

A. 7 

Damage zone 
model 

Uses Dugdale approach to represent damage; for pre- 
dicting strength of holes, cutouts, and cracks. 

Backlund, 

1981 

A. 8 

II 

Finite element 
method 

Uses 3-D finite element to determine interlaminar 
stresses; for delamination growth prediction. 

Ratwani, et 
al., 1979 

A. 9 

1-D thin-film model 

Uses laminate buckling theory for through-the-thickness 
delamination prediction. 

Chai et al., 
1981 

A. 10 

Edge-delamination 

model 

Uses rule-of-mixture for onset of edge delamination 
prediction. 

O’Brien, 1982 

A.1 1 

2-D thin-film model 

Uses laminate buckling theory for embedded delamina- 
tion prediction. 

Shyprykevich 
etal., 1984 

A. 12 


The analytical methods in each category of Table A-l have been evaluated. The advantages and 
disadvantages of each approach are summarized in Table A-2. In general, none of the existing 
models provides reliable results at reasonable cost. This is particularly true for category II 
models. The finite element approach is ideally suited for analyzing composite structures that 
contain delaminations, but computer costs are often prohibitive. 

The following selection criteria were established to distinguish which of the available models 
should form the basis of the required methodology: 

Models shall be generic in nature. 

• Lengthy calculations shall be avoided. 

• Models shall provide predictions with adequate accuracy and conservatism for space 
hardware application. 

• Models shall use the existing composite material database. 

• Methodology shall be readily adaptable to existing Aerospace computer systems. 

On the basis of these criteria, the IFM developed by Waddoups et al. (Ref. A.2) was selected for 
use in predicting the residual strength of cracked composite laminates. The PSM developed by 
Whitney and Nuismer (Ref. A.3) was selected for use with laminates containing holes and 
cutouts. These models require simple calculations and have the advantage of an extensive test 
database for the more common material systems such as fiberglass/epoxy and Gr/Ep. 

For composites with through- the- width or embedded delaminations are characterized by the thin 
film model of Chai et a/. (Ref. A. 10). For delaminations at free edges, the edge delamination 
model of O’Brien (Ref. A. 1 1) was selected. These two models have two common features: they 
avoid lengthy finite element calculations, and they use the strain-energy release rate G to 
characterize delamination growth. 
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The previously mentioned models form the basis for an Aerospace-developed damage tolerance 
software package called ADAM-C, which was the primary product of the Aerospace Sponsored 
R program. ADAM-C architecture is illustrated in Figure A-l. Material properties, laminate, and 
flaw geometries are input into preprocessor modules (GENLAM and RANKHO), which generate 
the laminate constitutive relations and the stress concentration factors for certain flaw shapes. 
These codes were acquired from the Air Force Materials Laboratory and linked to the remaining 
modules such that unnecessary programming effort was avoided. 



Figure A-l. ADAM-C methodology. 


Once the stiffness matrices and stress concentration factors have been calculated, control is 
transferred to the appropriate module, depending on the type of flaw selected. The program 
addresses holes, cracks, and elliptical cutouts. In each case, a residual strength prediction is the 
final product. The safe-life and residual buckling strength of a composite laminate that has a 
delamination may also be obtained by running ADAM-C. 


ADAM-C is available on both the VAX and IBM-PC systems. Each model is independent and 
modular, thus facilitating any future enhancements. The program is menu-driven and user- 
friendly. The material database contains many commonly used graphite and glass/epoxy systems, 
but this list can be updated by users at any time by simply running an interactive program called 
MATRPR (Ref. A. 13). Thus, any new material systems can be quickly accommodated. 


The following list summarizes the current ADAM-C output: 


Laminated plate theory stiffness matrices 

• Stress concentration factors due to holes or cutouts at arbitrary angles of inclination; 

• Static residual strength of cracked panels; 

• Static residual strength of panels with holes and cutouts; 

• 1-D and 2-D delamination growth with residual buckling strength; 

• Edge delamination growth and residual strength. 
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Table A-2. Assessment of Analytical Methods 


Loading and Damage 
lyee 


Analytical Method 


Advantages 


Disadvantages 





Point stress failure 
model 


Notched fracture stress 
model 

Damage zone model 


Inherent flaw model 

Strain energy density 
factor method 

Finite element approach 
Damage zone model 


Finite element method 


Edge delamination 
model 

Finite element method 


Thin-film 1-D model 


Finite element method 


Thin-film 2-D 
delamination model 


Simple calculation, ade- 
quate database 


Simple calculation 


Provides more accurate 
predictions 


Simple calculation, ade- 
quate database 


Provides more accurate 
predictions 

Provides adequate accuracy 
for complex geometries 

Provides more accurate 
predictions 


Provides adequate accuracy 
for complex geometries 

Simple calculation 


Provides adequate accuracy 
for complex geometries 

Simple calculation 


Provides adequate accuracy 
for complex geometries 

Simple calculation 


Stress intensity zone size is 
independent of hole sizes 

Extensive testing, limited 
database 

Length calculation, not enough 
database 

Characteristic zone size is inde- 
pendent of crack sizes 

Very limited database 

Length calculation 

Length calculation 


Length calculation involved 
Limited application 

Length calculation involved 
Limited application 

Length calculation involved 

Prediction accuracy needs to be 
improved 


Certain refinements were incorporated into the selected ADAM-C modules. In the PSM, for 
example, a finite width correction factor that accounts for the interaction between the hole and 
the free edge of the plate was added. The importance of the correction function is illustrated in 
Figure A-2. Experimental residual strengths for a Gr/Ep laminate are compared with theoretical 
values obtained for different combinations of damage zone sizes and finite width correction 
factors. It is apparent that only the curve that represents the combined effect of the damage zone 
and the finite width correction factor correlates most closely with the experimental data. Thus, 
the incorporation of an appropriate correction function is vital to the residual strength prediction 
of the model. 
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Figure A-2. Point stress model parametric study. From Reference A.8. 


The IFM is based on the linear elastic fracture mechanics model for metals containing cracks, but 
was modified to include the damage zone size. This expression depends on such parameters as 
the crack length, damage zone size, unnotched laminate tensile strength, finite width correction 
function, and most importantly, the laminate fracture toughness K q . 

Unlike metals, there are no standards for the definition and experimental determination of K q for 
composite materials, i.e., the fracture toughness is open to interpretation. Therefore, the 
calculation of K q by analytical techniques is an important topic, as we have already established 
that experimental determination of the important parameters of all candidate laminates in a given 
structure is not economically feasible. At the present time, our methodology uses two approaches 
to calculate the fracture toughness because no single model adequately covers the wide range of 
laminate behavior. The first approach assumes a simple relationship between the fracture 
toughness and the unnotched tensile property of a laminate, i.e., K q , is one-half of the unnotched 
strength (in English units). The second approach uses an analytical prediction of K q proposed by 
Poe (Ref. A. 14), who attempted to develop a general fracture toughness parameter that depends 
on the lamination constants as well as certain assumptions about the importance of load carrying 
plies. For conciseness, the first K q approximation is abbreviated as IFM-A and the second as 
IFM-B. 


Results of the delamination studies indicate that the strain-energy release rate G is a preferred 
parameter from the standpoint of delamination growth. An approach based on G avoids the 
necessity and expense of calculating stress distributions near flaw boundaries. Often, these 
stresses become singular and thus of no practical use in predicting stable delamination growth. 

Delamination growth requires work that may be supplied by external forces during delamination 
or by the decrease in the strain energy of the system or both. For thin film models, energy is 
supplied by the change in strain energy resulting from delamination extension. If the energy loss 
is greater than the work required to create a unit area of delaminated composite, growth will 
occur. The delamination growth is assumed to follow a power law based on the strain-energy 
release rate raised to some exponent. This approach has been incorporated into ADAM-C. 
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In order to substantiate the predictive capabilities of the ADAM-C program, an extensive 
correlation study was performed. Test data found in the open literature were correlated with 
predictions by means of the PSM and IFM modules. Although appreciable data are available on 
laminates with holes, cutouts, and cracks, information about delamination growth under cyclic 
loads is very sparse. Our future experimental work will seek to fill this data gap, and correlations 
will be made with the ADAM-C delamination growth modules. 

In Figures A-3 and A-4, the predictive accuracy of the PSM and IFM is correlated with some 
typical experimental data found in the literature. All of the data, both holes and cracks, are 
compared with three models: PSM, IFM- A, and IFM-B. The purpose of this approach is to 
investigate the extent to which the hole and crack models are interchangeable and thus develop 
some insight into the nature of these flaws in composite structures. Briefly, the results indicate 
that, for cracked structures, the IFM-B approach is preferred, whereas for composites with holes, 
the PSM provides the desired conservatism. Laminated structures have such a wide range of 
behavior that no single model can be completely endorsed, but it is believed that these models 
have the potential to be extended to more general types of loading conditions. 



The correlation between the PSM and typical data for laminates with holes is shown in 
Figure A-3. The selected database represents a wide range of laminate configurations and 
sources within the open literature. It should be noted that both axes have been normalized with 
respect to unflawed laminate strengths such that all data can be plotted as percentages of unity. 
Perfect correlation between theory and experiment is achieved along the darkened line with a 
slope of 1. The data points below the curve indicated conservative predictions. In general, the 
PSM provides adequately conservative predictions over most hole sizes of practical interest. 

In Figure A-4, IFM-A and B are correlated to some of the existing cracked laminate data. Again, 
the data have been normalized and the six symbols represent the three data sources and their 
comparisons with the two approaches for determining K. The scatter in the residual strength 
predictions is greater for cracks than for holes. 
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Material: T300/520B Graphite /Epoxy 
Laminate: [90/QJ2S 
Width: 3.9 in. 



Normalized Crack Size. 2a /w 


Figure A-4. Crack data correlation. From Reference A. 13. 

Further insight into each model is provided in Figures A-5 and A-6, where frequency of 
occurrence is plotted versus the ratio of the predicted-to-experimental residual strength values for 
the PSM and IFM calculations. A perfect model would always duplicate an experimental result, 
and hence the strength ratio (X-axis) would always be 1 ; there would be no distribution. 
Theoretical models are seldom perfect, however, and a reasonable basis for comparison is to 
cross-plot all distributions on the same figure. A good model should peak near X = 1 and should 
have a low-frequency response in the unconservative regime, i.e., X > 1. 



Experimental Normalized Residaal Strength, F/Ftp 


Figure A-5. Point stress model correlation for plates with holes. 
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Figure A-6. IFM correlation for plates with cracks. The data points below the curve indicate conservative 

prediction. 

In Figure A-7, the peak PSM response is in the range of 0.9 to 1.1. The IFM models also perform 
well, but are more frequently unconservative than is the PSM. Therefore, these findings indicate 
that the PSM is the preferred ADAM-C module for analyzing laminated structures containing 
holes. 



Residual Strength Ratio, pre die led /experimental 

Figure A-7. Frequency correlation (hole data). Less than 1.0 indicates conservative prediction. 

In Figure A-8, the three models are cross-plotted with respect to the cracked laminate data. The 
PSM is slightly unconservative. IFM-A yields a relatively flat response with both very 
conservative and unconservative answers. IFM-B results are sometimes too conservative but 
seldom unconservative. The broad scatter indicated in Figure A-8 suggests that it is difficult for a 
single approach to cover the residual strength variations of cracked laminate data. 
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Figure A-8. Frequency correlation (crack data). Less than 1.0 indicates conservative prediction. 

In summary, this project addresses the need for a reliable and cost-effective damage tolerance 
methodology for composite spacecraft structures. During the project, a literature survey was 
conducted, model selection and enhancements were made, software was developed, and 
comparisons were made between experimental and analytical results. These tasks resulted in the 
development and implementation of the ADAM-C code. This work is the first step in an effort to 
characterize the damage tolerance capabilities of composite structures in space environment. 

Now that the basic methodology is in place, the process of final refinement begins. As stated 
earlier, the aircraft and aerospace structural requirements can be quite different. Generally 
speaking, aircraft structures operate in a low-stress, high-cycle environment, whereas the reverse 
is often true for spacecraft components. The latter condition is more severe, but the existing test 
data were generated for aircraft applications in the low-cycle, high-stress regime. Therefore, an 
accurate correlation and refinement of ADAM-C will require further experimental testing under 
conditions relevant to spacecraft structures. 

A damage tolerance testing program is planned with specimens to be constructed from 
commonly used composite material systems. Flaws and defects will be prefabricated in the 
specimens to represent a variety of damage conditions. The results of these tests will be 
correlated with ADAM-C, and appropriate refinements to the methodology will be incorporated. 
With this approach, the goal of developing reliable damage tolerance analysis techniques with 
respect to spacecraft structures will be achieved. 
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Appendix B: PFA Code for Composites 


The research and development for the PFA of composites at Aerospace has integrated and 
extended the PFA developed at NASA LaRC in the late 1990s and early 2000s (Ref. B-l). The 
PFA is implemented in a software tool that integrates with existing commercial FEM software. 
The tool simulates the growth and propagation of failure modes in a composite and it is 
constantly adapted for new classes of problems and extended to incorporate new physics. The 
key features of the PFA code include: (1) Delamination simulation by using interface elements; 
(2) Prediction of in-plane and out-of-plane nonlinear shear damage; (3) Prediction of matrix- 
cracks based on a modified version of Puck’s criterion (Ref. B-2, B-3); and (4) A fiber failure 
criterion based on the work of Chang (Ref. B-4). The PFA code utilizes concepts from damage 
mechanics to degrade the material properties of the composite appropriately. The code also 
handles contact in the event of delamination closure following growth. Similar to ADAM-C, the 
code was developed to be as generic as possible and applicable to a large range of problems with 
adequate accuracy and conservatism. The accuracy of the code, however, relies on developing a 
complete mechanical, strength, and fracture database that can be used directly as input to the 
code. The next subsections will discuss the progressive failure philosophy, the damage models in 
the PFA code, the characterization testing that is required to define the damage model, and 
validation of the PFA code. The PFA code is also demonstrated for few examples pertaining 
pressure vessels. 

B-l Progressive Failure Philosophy 

Refer to Figure B-l for the following discussion on progressive damage philosophy and its 
relation to validation. The goal is to develop a reliable and practical analytical tool capable of 
simulating the response of composite joints and predicting their structural margins. This effort 
can be broadly divided into three sections: coupon testing, component testing and analysis. 
Coupon testing is essential in the damage modeling of composites since it is used to determine 
material properties, fracture properties and failure criteria. Coupon tests can also give insight into 
the progression of damage and response of the structure based upon observations and 
measurements. From component testing, the second leg in the model development, the structure 
can be evaluated on a macro-scale to evaluate margins. 

However, insight into the progression of failure is more difficult since failure may be unstable 
and the testing cost is considerably higher for flight components than it is for coupons. Finally, 
analysis is conducted to validate the material model. Since the coupon tests are simpler than the 
flight component, it follows that the analytical model will first be used to predict the response of 
the coupon tests. A FEM is created to reflect the experimental set-up; material properties are 
carried over from coupon testing or found in the literature, as appropriate, and care is taken to 
ensure that the loads and boundary conditions reflect the actual test configuration. The next step 
in the process is to compare predictions with test observations - including the measured strain, 
delamination formation and growth, shear damage, matrix cracking or fiber failure - in order to 
correlate the relevant material parameters. This is especially important in the modeling of 
damage, whereby friction coefficients and even strength the finite element results and the test 
data, the model can be used to analyze a subscale structural component and compared to 
experimental results. If large differences exist between the finite element results and the test data, 
then the model needs to be re-examined. Typically, the material, boundary conditions, and 
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failure criteria would be revised at the coupon level. This refinement continues until the subscale 
analysis shows good agreement with the subscale test data. The analytical tool is considered to 
be reliable once the model can accurately predict the response of both the coupon and component 
tests. 

Figure B-l illustrates the modeling technique and adjust for specific applications. Both coupon 
and component testing are required, along with corresponding finite element analyses. 



Figure B-l . Modeling flowchart. 


B-2 Progressive Failure Methodology 

A progressive failure methodology is developed that is an extension of the methodologies found 
in Pifko and Kushner (Ref. B-5) and Sleight et al. (Ref. B-6]). The methodology requires a 
nonlinear solver to establish equilibrium because both geometric and material nonlinearity are 
considered. The geometric nonlinearity is due to large strain and large rotation kinematics. The 
nonlinear material behavior is due to the degradation of the material properties of the plies and 
the material in between the plies to simulate intralaminar and interlaminar damage mechanisms. 
An incremental iterative approach is adopted for the nonlinear finite element analysis, and the 
Newton-Raphson method is used to trace the loading path of the structure with a displacement- 
control analysis. A line search algorithm is only activated to provide an initial approximation to 
the Newton-Raphson method, or to minimize a large residual force. The progressive failure 
methodology complies with two requirements: (i) optimum size load steps to achieve solution 
convergence while representing accurately local load redistribution, and (ii) a fine mesh in 
regions that undergo material degradation to represent proper stress gradients at the boundary of 
the damaged zones. 

Damage parameters used in the analysis can be divided in two parts: one relates to the failure 
prediction at the lamina level, and second relates to the prediction of damage initiation at the 
laminate level, which can then progress and lead to ultimate failure. The progressive damage 
implementation is depicted in Figure B-2. The damage evolution equations are solved at the 
Gaussian integration points within each material ply, thus explicitly accounting for damage 
variation through the thickness and in the plane of each material ply. The flowchart is color- 
coded as follows: Green represents the stage at which the user interfaces directly with the finite 
element software ABAQUS™ (Ref. B-7); Blue represents the user material subroutine UMAT; 
and Peach represents the ABAQUS™ standard nonlinear incremental procedure. First, the 
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ABAQUS™ input file needs to be prepared in a typical fashion with additional parameters 
relating to the strength and fracture of the composite. Then, the analysis is initiated by applying 
an initial load to the FEM. A converged solution must be obtained because the structure is no 
longer in equilibrium due to differences between the internal loads and the external loads. The 
stress, strain, and stiffness matrix (Jacobian) are updated. At each load step, failure criteria are 
applied at each Gaussian point. If failure occurs, the damage state is obtained and the 
components of the stiffness matrix are reduced correspondingly. The components are reduced to 
a small value instead of zero, because a reduction to zero leads to numerical difficulties in the 
nonlinear procedure. Following the application of the damage tensor, the modified Newton- 
Raphson nonlinear iterative procedure is conducted to achieve equilibrium prior to the next load 
step. 



Figure B-2. Flowchart of the progressive damage implementation using ABAQUS™. 


The progressive failure methodology was implemented in a user subroutine package compatible 
with ABAQUS™ finite element code, Figure B-3, because it has the ability to solve large 
models, has advanced material models and contact definitions, and it has a built-in post 
processing software. 


Delamination 

model 


Fiber failure 


User Subroutine 
Package 


Matrix-cracking 

model 


Nonlinear 
shear model 


Ability to model over 
500K DOF 


Efficient nonlinear solvers 
with parallel processing 
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Stress recovery from 
discrete composite layers 


Damage visualization using variety 
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ABAQUS 
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state of the art in structural analysis 


Model Subsystems Integration Benefits 

Figure B.3. Failure modes integrated into a user subroutine package. 
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B-3 Failure Modes and Degradation of Mechanical Properties in PFA 

A brief description of following failure modes considered in the progressive failure methodology 
is presented: (1) Delamination; (2) In-plane and out-of-plane nonlinear shear damage; 

(3) Matrix-cracking; and (4) Fiber failure. 

Delamination 

A new computational fracture mechanics technology has been developed in recent years 
(Refs. B-8 and B-9) for the prediction of discrete delaminations in composite materials via the 
use of interface elements coupled with cohesive zone models. The interface element formulation 
consists of a nonlinear distribution of continuous springs, Figure B-4, with a cohesive- 
decohesion constitutive law that satisfies a strength based failure criterion and a fracture-based 
criterion. The interfacial constitutive law for the formulation of the interface element governs the 
deformations of these interface layers and determines whether delamination growth occurs or 
not. This unified strength-fracture approach allows the prediction of the initiation and growth of 
delamination. Interface elements are positioned at interfaces where the delamination may 
potentially grow. The prediction of multiple delaminations and non-self similar delamination 
growth is then obtained depending upon the natural response of the structure being analyzed. 



Deformed 

Figure B-4. Interface element consists of a continuous distribution of breakable springs. 

Prior versions to ABAQUS™ 6.5-6 did not have the interface finite element formulation to 
predict delamination, and thus the element was previously developed and implemented in 
ABAQUS™ using a user element subroutine (UEL). The drawback with the UEL approach is 
that ABAQUS™ does not support visualization of the damage distribution at the interface 
element level, so the shape and size of the delamination is not known, except for separations that 
are visible from the deformation plot. To overcome this visualization difficulty, an effort was 
undertaken to predict delamination using standard, but very thin, solid brick elements (C3D8I) 
available in the ABAQUS™ element library. The traction-separation law normally 
accompanying interface elements was transformed to a stress-strain law for an interphasic 
element of very small thickness and implemented in ABAQUS™ with a user material subroutine 
UMAT. Severe convergence issues surfaced occasionally, mainly triggered by excessive element 
distortion. The element distortion problem occurred more often under the initiation and growth 
of shear delamination. Moreover, the formulation of the C3D8I element was not meant to be 
used for the prediction of delamination. The most current version of ABAQUS™ has added the 
capability to predict delamination using interface elements. 
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The interface element type used in this paper is the COH3D8, an eight-noded element, three 
displacement degrees of freedom per node, and linear interpolation. There are two convergence 
issues with the COH3D8. The first issue is that the convergence has been found to be excessively 
slow with the built-in bilinear constitutive law that accompanies this element formulation. The 
source of the slow convergence is the tangent stiffness discontinuity present in the bilinear 
constitutive law. The Newton-Raphson nonlinear solver enters into nonconvergent oscillating 
cycles of the residual force while attempting to converge to a solution. The second issue for slow 
convergence is related to numerical instabilities arising from stiffness matrices with negative 
eigenvalues and ill-conditioned stiffness matrices that often accompany the type of problems 
involving softening behavior. 

The first convergence issue is addressed by selecting a constitutive law with a smooth slope, 
Figure B-5. The constitutive law has an exponential form, and it naturally satisfies the following 
multi-axial quadratic stress criterion for the onset of delamination: 


f rri \ 

l 

r 7j ^ 

2 

f rr, \ 

1 2 


+ 

+ 


UJ 




U ) 


( rp \ 2 

±3 

v^3c / 


+ 


V 5 i 


t/T 3 


\ 2 f 

+ 


v • 


t/T 3 


= U < 0 ; 


where T\ and T 2 are the interfacial traction components associated with shear and Si and S 2 are 
the interfacial shear strengths. T 3 is the normal interfacial traction component, S 3t is the tensile 
interfacial normal strength, and S 3c is the compressive interfacial normal strength. The quadratic 
failure criterion includes a ‘friction coefficient’ that incorporates the apparent increase in shear 
strength due to transverse compression. The constitutive law also naturally satisfies the following 
mixed-mode fracture criterion for the progression of delamination: 

G \ f % f fin _ j. 

G\c Cue Gm 


where the Mode I, II, and III energy release rates are Gi, Gw, and Gm, respectively; and the 
critical Mode I, II, and III energy release rates are G\ c , Gu c , and Gnic, respectively. Naturally 
means that the criteria do not need to be evaluated in the numerical algorithm, because the 
criteria and its corresponding failure degradation rule is built-in in the constitutive law. These 
constitutive equations are also thermo-mechanically consistent, and this is achieved by including 
a damage parameter to prevent the restoration of the previous cohesive state between the 
interfacial surfaces. The second convergence issue is resolved by setting negative diagonal terms 
in the material tangent stiffness matrix to zero. This modification does not affect the solution, but 
it improves convergence efficiency. The constitutive law implementation is conducted with a 
user material subroutine UMAT. The constitutive law requires definition of the critical energy 
release rate for Mode I, II, and III; and the shear and normal interlaminar strengths. 
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Figure B-5. Interface element coupled with a cohesive zone model. 
In-Plane and Out-of-Plane Nonlinear Shear 


Composite laminates typically exhibit a nonlinear response when subjected to shear loads either 
in the plane of the composite or transverse to the plane of a thick composite. Material 
nonlinearity in the plane of the composite can be replicated during off-axis testing or Iosipescu 
testing. A typical nonlinear curve is shown in Figure B-6. Testing is generally used to extract the 
initial shear modulus and the shear strength. Hahn and Tsai (Ref. B-10) have proposed a cubic 
polynomial to represent the nonlinear in-plane shear stress-strain relationship. 



Figure B-6. Typical nonlinear shear constitutive law. 

The work at Aerospace has been extended to account for nonlinear shear through-the -thickness 
of the composite. The extension was conducted because there is extensive evidence presented in 
the literature indicating that the nonlinear shear model plays an important role in accurate 
predictions of the failure strength. Including this nonlinearity is of paramount importance to find 
the appropriate stress distribution and concentration at the bolt hole. Otherwise, premature 
predictions of in-plane matrix-cracks are predicted and may result on a gross underprediction of 
strength. Including the nonlinear shear behavior of the composite in the failure modeling relieves 
the in-plane shear stress concentration at the hole boundary and delays the prediction of matrix- 
cracks. Although complex polynomial forms and sophisticated nonlinear models have been 
proposed, at Aerospace the following exponential constitutive law was developed that requires a 
minimum number of material parameters: 
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The material properties that need to be specified for the nonlinear shear constitutive law are the 
initial shear modules, ultimate shear strengths, ultimate shear strains, and the yield strengths. 

Matrix-cracking 

The strength-based failure criteria developed by Hashin (Ref. B-l 1 ) based on experimental 
observations of unidirectional composites subjected to tensile loads, has been extended because 
it does not always fit experimental data as shown by numerous studies, especially in the case of 
matrix compression. The World-Wide Failure Exercise (WWFE) was conducted, whereby the 
leading theories for predicting failure in composite laminates were compared against 
experimental evidence in order to assess the qualitative and quantitative robustness and accuracy 
of each failure theory. The comparisons included biaxial strength envelopes for a range of 
unidirectional and multi-directional composites. The WWFE results indicated that even for 
simple layups and loadings, many of the failure criteria that have been proposed to-date fall short 
of predicting the failure strength accurately. Since the numerical predictions with the failure 
theory by Puck and Schurmann (Ref. B- 3 ) showed good agreement with experimental results for 
the prediction of biaxial failure envelopes in the WWFE studies, their criterion is extended from 
in-plane only to a 3 -D form that accounts for through-the-thickness failure. The failure theory 
follows the work of Hashin. Hashin’ s failure criteria for matrix-cracks consists of a quadratic 
interaction of stress invariants, whereby the criteria is made independent of the stress parallel to 
the fiber based on the fact that any possible failure plane is parallel to the fibers. The failure 
criteria proposed by Puck and Schurmann is essentially an extension of Hashin’ s postulate. 

The fracture plane is unknown and it must be determined. This is accomplished by maximizing 
the failure index with respect to the fracture plane angle. Since matrix-cracking occurs in a plane 
parallel to the fibers, the fracture plane oriented at an angle /?with respect to the global 
coordinate system must be obtained. The tractions acting in a potential fracture plane are 
determined by force equilibrium as follows: 

<7 nn = — + 2 — + — 2 — cos(2/7) + ct 23 sin(2y0) 

G nt = - 1722 - a33 sin(2 /?)+ <T 23 cos(2/?) 

<7nl = CT \2 cos(/?) + <7 13 sin(yff) 

where Oijis the local state of stress that rotates with the lamina orientation. < 7 n is the stress 
parallel to the fiber, <733 is the stress in the stacking direction, On is the in-plane shear stress, and 
Ob, <7)3 are the transverse shear stresses. The matrix-cracking failure criterion for a tensile 
traction is postulated as 
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and for a compressive traction as follows: 
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The contribution of the compression traction has been eliminated but added to increase the shear 
strength in order to indicate that shear fracture is impeded in compression. In compression there 
is additional resistance to shear fracture. This effect is incorporated with a ‘friction coefficient’. 
Here, the denominators in this failure index can be regarded as effective shear strengths. The 
potential fracture plane as shown in Figure B-7 is obtained by maximizing the failure index. 
Once failure has been established, corresponding mechanical properties are degraded to emulate 
fracture. 



Figure B-7. Potential fracture plane in the 3-D failure criterion for matrix-cracking. 


Fiber-breakage 


The failure criteria for the prediction of fiber failure utilizes critical strain and it is simply: 

e n < c £n, £n<0;£n> Vn, £ii > 0 

where the c £\ \ is the strain allowable in compression and l £\ \ is the strain allowable in tension. 

For fibers surrounding a hole, an area satisfying the failure criterion is used as the criteria. The 
need for a critical fiber area criterion in the case of a hole is made evident by stress 
concentrations that significantly exceed the fiber strain allowable, thus causing rapid damage 
growth at load levels that are not realistic. Utilizing the critical fiber area as a criterion delays the 
failure of the composite with a notch or hole. Otherwise, the strength could be under predicted by 
more than 50 percent. 


Degradation to Mechanical Properties 


Now that the failure mechanisms with their respective failure criteria have been defined, the 
constitutive law for the composite material needs to be modified accordingly to simulate damage 
initiation and growth. The new constitutive law, as derived from damage mechanics is: 


C ij = C ij (l-d i )(l-d j ) 


where i = 1 ,2,3 correspond to the direction of the fiber, to the direction perpendicular to the fiber 
in the plane of the laminate, to the direction perpendicular to the fiber transverse to the laminate, 
respectively. 

Fiber breakage corresponds to di = 1 and the components of the damage elasticity tensor Cij,Cji 
become zero. The physical interpretation of the reduction in stiffness C 13 , C 31 , C 21 , C 12 is that if 
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the fiber fails, there is no load transfer from the fiber to the matrix. From the engineering 
standpoint, this reduction represents that the Poisson ratio no longer affects the mechanical 
behavior of the composite. Similarly, for matrix cracking in the plane of the laminate and 
transverse to the laminate d 2 = 1, and the components of the damage elasticity tensor C 2 j,Cj 2 
become zero. If matrix-cracking occurs, the shear stiffness for each shear plane is also degraded 
using the damage parameter d 2 . For delamination, the energy dissipated by this failure 
mechanism will be taken into account using interface elements, hence the damage variable d 3 
will be held as zero. 

B-4 Material Characterization for the PFA 

To properly simulate damage initiation and growth of different failure mechanisms and their 
interaction, the mechanical, strength, and fracture properties must be determined from coupon 
testing. Table B-l shows basis material properties determined from coupon testing. The type of 
tests needed are simple tensile test, Iosipescu testing, and the open-hole strip testing. For fiber 
failure and matrix-cracking coupon testing can be used to fully characterize failure in large scale 
structures, Figures B-8 and B-9. From the tensile test, one can determine the Young modulus and 
the ultimate tensile strength in the fiber direction. The shear stiffness, the yield stress, and the 
ultimate shear strength in the plane of the laminate can be determined either from [— 45/+45] ns 
off-axis loading test, short-beam shear test, or torsion test. From the open-hole strip testing, one 
can determine the critical area for which neighboring fibers fail. 

The fracture properties to characterize delamination growth are obtained from fracture testing 
that are either an American Society for Testing and Materials (ASTM) standard or are in the 
process of being certified. The Mode I critical energy release rate is determined from the double 
cantilever beam. End loads are applied in opposite directions to peel the composite beam in to 
two sublaminates. The test is generally performed on unidirectional composites with the fibers 
oriented such that these are parallel to the length of the initial delamination. The experiment 
consists of a load that is applied at aluminum end blocks jointed to the DCB specimen by an 
adhesive. The test is performed under displacement-controlled tensile test machine. The 
conventional quasi-static tests are generally performed at constant displacement rates that can be 
chosen between 5 mm/min and 15 mm/min. 

Table B-l. Basic Material Properties Determined from Coupon Testing 
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Figure B-8. Test configurations for the characterization of fracture properties. 



losipescu Test Modified Arcan Test 



DERA 



Figure B-9. Test configurations for the characterization of strength properties. 


The initial delamination is achieved by incorporating a Teflon film at the midplane of the 
laminate layup prior to curing. There are issues pertaining the test specimen layup. With these 
test specimens, the delamination front may not be straight due to uneven distribution of the 
energy release rate across the width of the beam. This uneven distribution is caused by the 
anticlastic bending effect. Secondly, the main delamination may branch into multiple cracks that 
may follow the fiber-matrix interfaces. This situation may lead to Mode II loading. In addition, 
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the phenomenon of fiber bridging caused by fiber nesting in 0-degree composites can lead to 
variations of the critical energy release rate. 

Table B-2 summaries fracture test configurations. A complex fracture behavior involving fiber 
breakage, ply jumping, and fiber bridging often occurs when arbitrary orientation in ply angles 
are considered in the design of DCB test specimen. The delamination branches to interfaces 
away from the midplane, which leads to larger values of fracture toughness. They developed a 
DCB test method to suppress crack jumping and fiber bridging effects. To achieve pure Mode I 
delamination in a layup other than 0-degrees, the arms of the DCB specimen were designed such 
that these are balanced and symmetric to eliminate the stretching-shearing and stretching 
coupling effects. In addition, the layup of the angle ply laminate was designed to minimize the 
bending-twisting effects. The tests were designed to ensure there was no curvature or shear 
distortion of the laminate due to thermal stresses from curing. 

The critical energy release rate for Mode II is determined either from the ENF, end load split 
(ELS), mixed-mode bending (MMB) apparatus, or the four-point ENF. Only the ENF and the 
ELS will be discussed because these are the most commonly used in determining the fracture 
energies under Mode II loading. The ENF is simply supported at the ends with the load applied 
at the middle of the specimen. If the delamination length is approximately less than 35 percent of 
the total length, the test specimen undergoes delamination. 

The ELS test configuration is clamped at one end and a load is applied at the other end. For this 
test configuration, the delamination growth is unstable for short delamination lengths that are 
less than 55 percent of the length of the specimen. In both test configurations, Mode II loading is 
promoted by the relative sliding of the upper surface of the lower lamina with respect to the 
lower surface of the upper lamina. The relative sliding occurs due to the rotation of the upper and 
lower arms. The practical difficulties in measuring the Mode II fracture toughness with these test 
configurations is the choice of starter defect, the stability of the test, and frictional effects 
between the crack faces. Using a film insert as a starter defect leads to a larger fracture toughness 
than a test specimen with a pre-crack. The main problem encountered with the Mode II tests is 
that determining the stability of short cracks is complex since these can dynamically propagate. 

The critical energy release rate for Mode III can be determined from the end crack torsion 
(ECT), the split cantilever beam, and the crack rail shear. The split cantilever beam generates a 
substantial Mode II strain-energy release rate component near the free edges, thus invalidating 
the test. Extensive efforts undertaken by NASA LaRC to standardize the ECT test to characterize 
Mode III have been unsuccessful. The edge-crack torsion test configuration consists of a 
rectangular composite of a [90/(±45) n /(±45) n /90] s laminate with a rectangular pre-crack 
introduced by a non-adhesive film inserted at the midplane along one edge during manufacture. 
The sample is mounted in a special fixture that provides support at three points, and a transverse 
load is introduced at the fourth point by compressing a loading pin. Usually, the specimen and 
fixture assembly is loaded in a test machine with a constant stroke rate of 1.5 mm/min. This 
special fixture creates torsion along the length of the laminate, thus inducing Mode III 
delamination. 

The MMB test, crack lap shear, or fixed ratio mixed mode (FRMM) can be used to determine the 
failure envelope for a delamination under a mixed-mode condition. From all the tests discussed, 
the double cantilever beam and the MMB tests are ASTM standards, while the end notch flexure 
is a European Agency standard. Extensive efforts are underway to develop the other fracture test 
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configurations into ASTM standards fully. For the full characterization of the transverse strength 
properties including the nonlinearities required in the PFA, the Iosipescu test, Modified Arcan 
Test, and the off-axis compression tests can be used to determine the interlaminar shear strength, 
tensile strength, and the mixed-mode interaction. The benefit of the off-axis compression test is 
that it provides the means to determine the ‘friction coefficient’ required for the formulation of 
the interface elements. Typically, this ‘friction coefficient’ is 0.28 and it can be confirmed by 
testing. This ‘frictional parameter’ is essentially a phenomenological representation of physical 
events occurring at the microscopic level that are not well understood. Many researchers 
including Puck have treated this parameter as friction for matrix undergoing brittle fracture under 
compression. 


Table B-2. Summary Table for Fracture Test Configurations Analyzed Using Finite Elements 


Test Configuration 


Failure 

Modes 


Failure Description 


1. DCB 
Layups: [0] 2 


Mode I 


(1) Delamination front is curved. 

(2) For short cracks, delamination is dominated by 
strength, while for long cracks delamination is dominated 
by fracture energy. Unstable delamination growth occurs 
under load control. 


2a. ELS 


2b. ENF 


Mode II 


(1) Unstable delamination growth occurs either under 
displacement or load control if the crack length to 
specimen length ratio is 0.55 for the ELS and 0.35 for the 
ENF. 

(2) The energy dissipated due to friction is negligible 
compared to the fracture dissipation. However, frictional 
forces can cause temporary crack arrest. 

(3) For short cracks, a snapback algorithm is required to 

follow the load-displacement response. 


3.ECT 

Layups: [90/±45 3 / + 45 3 /90] s 


Mode III 


(1) Delamination is governed by (7 m , although a G\\ 
component does exist. Thus, delamination growth occurs 
non-self-similarly. 

(2) Proper determination of the transverse shear modulus 
(7 23 is critical for an accurate determination of the anti- 
plane shear mode. 

(4) Unstable delamination growth occurs if the initial 
delamination is sufficiently small. 


4a. FRMM 
Layups: [0] 2 


4b. MMB 


Mixed-mode 

failure 


a 


a 


(1) Delamination is governed by G\ and (7 n . 

(2) The delamination growth is unstable under either load 
or displacement control if the crack length to specimen 
length ratio is 0.41 for the FRMM and 0.35 for the MMB. 
Snapback required to trace the structural response. 

(3) The Riks algorithm required to follow the loading 
path of the MMB, since the loads are applied 
proportionally. (5) 
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B-5 Validation of PFA for Coupon Test Configurations 

To demonstrate the capability to predict delamination and the irreversibility capability of the 
PFA code, steady-state delamination growth of unidirectional T300/977-2 carbon fiber- 
reinforced epoxy laminate beams were simulated for quasi-static loading-unloading cycle of 
various fracture test configurations. The finite element results for Mode I double cantilever beam 
(Figure B-10), the Mode II ELS (Figures B-l 1), ENF (Figure B-12), and the mixed-mode tests 
(Figures B-l 3, Figure B-l 4) were in excellent agreement with either experimental data available 
in the literature or with linear elastic fracture mechanics analytical solutions. The PFA code was 
also validated with test data available for the ECT, a very unstable test configuration. The ECT 
test specimen is a G40-800/R6376 graphite epoxy composite laminate with a 
[90/(±45)3/(±45)3/90 ] s stacking sequence. The predicted failure load was within 5 percent of the 
test, Figure B-l 5, and the predicted damage distribution was in close agreement to the scans of 
the failed specimens, Figure B-l 6. Lastly, the finite element simulations of the DCB and the ELS 
made of AS4/PEEK were conducted and the numerical response was in close agreement with 
experiments, Figure B-l 7. 

The effect of through-the-thickness compressive/tensile loads on the strength and response of 
woven carbon/epoxy composite blocks loaded in an off-axis configuration with respect to the 
stack direction was investigated using the PFA code. This particular problem tests the ability of 
the code to consider the apparent increase in shear strength when a composite is subjected to 
through-the-thickness compressive stress. Published data suggests that the composite blocks fail 
primarily by delamination. Thus, to predict the strength of the blocks, decohesion elements were 
placed between adjacent lamina. This is an improvement over previous codes that relied on a 
simple quadratic strength-based failure criterion, thus, under-predicting the failure load. 
Compression simulations were conducted for 0-degree, 7.5-degree, 19-degree, 25-degree, and 
45-degree off-axis angles test configurations. Simulations were able to capture the catastrophic 
failure of the blocks, whereby the composite blocks crack into two pieces and the upper and the 
lower portions instantaneously rotate in opposite directions, Figure B-l 8. The failure event as 
observed from testing and predicted by analysis is dynamic and catastrophic. The predicted 
failure loads were in close agreement with the test data for all compression test configurations. 
The failure envelopes for the composite blocks were in close agreement with the testing 
conducted by Schubel and Sun’s failure criteria, Figure B-19. The PFA code also includes the 
nonlinear shear constitutive law, which enabled the prediction of the nonlinear structural 
response, Figure B-20. 
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Figure B-10. Predictions by PFA of the DCB compare well to experimental data and analytical solutions. 
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Figure B-l 1. Predictions by PFA for the ELS compare well to analytical solutions. 



Figure B-12. Predictions by PFA for the ENF compare well to analytical solutions. 
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Figure B-13. Predictions by PFA for the FRMM compare well to analytical solutions. 



Figure B-14. Predictions by PFA for MMB compare well to analytical solutions. 
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Figure B-15. Predictions by PFA for ECT compare well to testing. 
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Figure B-16. Damage predicted by PFA compared well by scans of the failed coupon. 
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Figure B-17. PFA predicted load-displacement response for DCB and ELS accurately. 
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Figure B-18. Predicted catastrophic failure of the off-axis test configurations. 



Figure B-19. Failure envelope predictions compared to experimental data for compression tests. 



Figure B-20. Failure response compared well with the experimental response. 
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B-6 Subscale test panels 

The PFA code is now demonstrated to simulate the initiation and growth of intralaminar and 
interlaminar failure modes occurring in AS4-3503 Gr/Ep panels with centrally located circular 
cutouts when subject to compression or shear loads, Ref. B-l. The testing of these panels was 
conducted at NASA LaRC. The failure of these panels is of special interest, because buckling, 
interlaminar, and intralaminar failures interact into the postbuckling regime of the structural 
response. The flat composite panel, Figure B-21, consisted of 24 plies with a stacking sequence 
of [±45/02]3 S and framed by a square fixture made of steel. Three comers of the frame were held 
by pin joints, while a displacement was specified in the remaining comer at 45 degrees with 
respect to the edge of the panel. 

Two curved composite panels one with diameter to width ratio of 0.4 and the other one with a 
ratio of 0.2, Figure B-22, were subject to axial compression. Both panels had a laminate stacking 
sequence of [±45/0/90]3 S . The PFA code was necessary to predict the structural response of these 
panels reasonably, Figures B-23 thm B-25. The extent and location of delaminations predicted 
by the PFA code was confirmed, Figure B-26, by visual inspection of the failed test specimens. 
The delaminations both as observed from test and analysis occurred at the midplane of the flat 
panel loaded in shear and the curved panel with d/W = 0.4 loaded in compression. For the curved 
panel with d/W = 0.2, delamination occurred in multiple locations throughout the thickness of 
the panel. Specifically, the PFA code predicted delamination initiation and growth perpendicular 
to the loading direction in regions where large compressive loads were present. As simulations 
identified the shear panel collapsing, a rapid progression of delamination perpendicular to the 
loading axis was predicted, while delamination arrest was predicted parallel to the loading axis. 



Figure B-2 1 . Flat panel loaded in shear. 
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Figure B-22. Curved panel loaded in compression. 
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Figure B-23. Structural response of the flap panel loaded in shear. 
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Figure B-24. Structural response of the curved panel loaded in compression for d/W = 0.4. 



82 


42000 


^>0 * C 

PFA, Present. Mid-surface 

PFA. Present . Bottom 457-45 

Experiment^ Hilbiirger el al_, 2001 

o PFA. Jaimky et al, 200 1 
■ Intralaminar failure initiation 
A DeUimiiiulinii milniticm 

— i ■ i 

0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 

End shortening, w (in.) 

Figure B-25. Structural response of the curved panel loaded in compression for d/W = 0.2. 


Figure B-26. Visual indications of delaminations correspond to delamination prediction locations. 

The PFA code also predicted delamination growth occurring perpendicular to the compressive 
loads for the two curved panels with different cutout sizes. Simulations indicated that the 
compression panel with smaller diameter cutout-to-width ratio failed with rapid progression of 
delamination from the edge of the cutout. For the compression panel with larger diameter cutout- 
to-width ratio, the predicted delamination growth was followed by delamination arrest. At this 
point, the numerical response indicates a residual strength, a finding consistent with the 
experimental data. After further loading, the predicted delamination growth leads to catastrophic 
failure. Thus, the PFA code was capable of predicting the rapid progression of delamination from 
the edge of the cutout at a location perpendicular to the loading axis including non-self-similar 
delamination growth, delamination arrest, and closure of delaminated faces. 

B-7 Impact Damage of Structures 

The PFA code was also used to determine the extent and shape of delaminations of a 
square sandwich panel with facesheets made from woven carbon fabric/epoxy laminates 
(AGP370-5H/3501-6S) and a closed-cell core made of PVC foam (Divinycell H250) subject to a 
low velocity impact caused by a free-falling blunt impactor, Figure B-27, Ref. B-12. The 
predicted delaminations were ringed shape, extended outwards from just the circular edge of the 
impactor, and were located in the midplane of the upper facesheet. The predictions compared 
well with ultrasonic scanning of the failed test specimen, Figure B-28. The predicted global 
response of the sandwich structure was in close agreement with the tests, Figure B-29. 
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Figure B-27. Sandwich structure impacted by a free-falling mass. 



Figure B-28. Strain response at bottom facesheet of the sandwich structure. 

Test 



Figure B-29. Delamination predicted by analysis compared well with the test. 


Two characteristics of the PFA code were instrumental in the good correlation that was achieved: 
it incorporates both the increase in interlaminar shear strength, Ref. B-13, and Mode II fracture 
toughness when the composite is subject to transverse compression loads. The simulation of this 
impact event demonstrated that considering other damage modes not included in the PFA code is 
important. A separate damage model for the foam core was necessary in order to properly 
simulate the extent and shape of delamination. The PFA code was also able to predict a decrease 
in load carrying capability caused by delamination growth, which was also observed in testing. 


The PFA code was validated with test data generated at NASA LaRC of a J-section circular 
fuselage frame made of textile perform 3M PR500 epoxy resin subject to radially inward load 
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representing a crash-type event, Figure B-30. The architecture of the textile composite is a 2^2, 
2-D triaxial braid with a yam layup of [0 degree 18k/±64 degrees 6k] 39.7 percent axial. The 
PFA of the frame correlated reasonably well with the predicted failure modes of the composite, 
including the structural response of the structure, Figure B-3 1 . The PFA code was numerically 
stable enough to predict the load-displacement response through several major failure events 
interacting with buckling. Since different failure modes, as predicted by analysis, occur 
simultaneously at different locations of the stiffener, the exact sequence from the model could 
not be extracted. However, in multiple locations kinking of axial tows led to matrix-crazing and 
eventually delamination. 



Figure B-30. J-section of a circular fuselage frame subject to a radial load. 



Displacement of Mid Point (in) 

Figure B-3 1 . Structural response of the J-beam subject to radial loads for the prediction and testing. 

B-8 Adhesively Bonded Joints 

The PFA code was validated against test data for several adhesively bonded joint configurations, 
Ref. B-14. The first validation was conducted with a single lap joint test configuration. The 
single lap joint test specimen, Figure B-32, was made of woven carbon epoxy fabric and epoxy- 
acrylate adhesive. The PFA code predicted two cracks initiating from the reentrant comers and 
propagated towards the middle of the test specimen. The PFA predictions compared well with 
the structural response from the experiments. 


85 






Figure B-32. Failure and response of the single lap joint test configuration. 

The second validation example was that of an Aluminum 6061-T6 DCB with a thermoplastic 
acrylic polymer adhesive, Figure B-33. Interface elements were placed within the adhesive to 
provide multiple paths for crack propagation and allow for the prediction of hackle crack 
patterns. The PFA code predicted crack paths similar to those observed with tests and compared 
well with the structural response from the experiment. Through the PFA it was found that a small 
component of Mode II fracture turned the crack and contributed to temporary crack arrests also 
observed in testing. 
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Figure B-33. Failure and response of a double cantilever beam test. 

The third validation example was that of a crack lap shear test, Figure B-34. The crack initiates 
in the vicinity of the reentrantcomer, or, at the material discontinuity as predicted by PFA and 
observed in testing. The PFA indicates that the crack is mainly dominated by a Mode II loading 
The predicted failure and response of the crack lap shear test configuration compared well with 
the testing, Figure B-35. 
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Figure B-35. Failure and response of a crack lap shear test. 


B-9 Mechanically Fastened Joints 


The PFA code was further validated for problems involving composite mechanically fastened 
joints, Ref. B-15 and Ref. B-16. The composites were 128 plies [0/-45/90/45]i6 S quasi-isotropic 
panels and 0.75 inches thick. The analysis simulated bypass only test, double-sided bearing and 
bypass load test, and single-sided bearing and bypass load test, Figure B-36. These tests were 
performed at Aerospace’s Space Materials Laboratory. The model incorporated nonlinear 
material behavior, and progressive damage of composites. 
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Figure B-36. Single-sided bearing and bypass load (top) and double sided bearing and bypass load 

(bottom). 


The orientation and density of matrix-cracks and fiber damage are represented through a damage 
tensor, an internal state tensor that is thermodynamically consistent within the framework of 
continuum damage mechanics. The damage tensor evolves as damage associated with different 
failure modes accumulates. As described earlier, matrix-cracks are predicted using a 3-D failure 
criteria, which is based on the stresses acting on a potential fracture plane. For use with the 
failure criteria, a new anisotropic shear constitutive law is postulated governing the transverse 
and the in-plane shear nonlinear mechanical behavior. 


The prediction of initiation and progression of multiple delaminations through the thickness of 
the composite is simulated via cohesive-decohesive constitutive law. Numerical simulations 
were in good agreement using filled hole tension, double sided bypass and bearing load, and 
single sided bypass and bearing test configurations. The predicted strength and the strain 
behavior was in good correlation with in-house experimental data, Figures B-37 thru B-39. 
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Figure B-37. Joint strength failure of various specimen configurations. 
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Figure B-38. Strain comparison for filled-hole tension test. 
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Figure B-39. Strain comparison for single-sided bearing testing. 
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B-10 Modeling Delaminations in COPVs 

Delaminations occasionally occur in COPVs. The objective of this section is to describe a 
method used to evaluate these defects by using a pre-processor in conjunction with progressive 
failure analysis. An example is used to illustrate this approach for a typical COPV. COPVs store 
highly pressurized gases in launch vehicles and spacecraft. Typical COPV construction consists 
of high-strength fiber composites wrapped around a thin metal liner. The purpose of the liner is 
to act as a leakage barrier while the composite carries the most of the load from pressurization. 
This combination leads to lighter weight designs for the same pressure capability, which is 
valuable for space applications. 

Structural integrity and conformance to safety requirements demand accurate stress, fatigue and 
fracture mechanics analyses. Typically, netting analysis has been used as a design tool for 
COPVs, but is a highly simplified method that cannot always be used to evaluate the structural 
integrity of a discrepant COPV. Accurate analysis is particularly important since the 
consequences of a local defect propagating can lead to catastrophic failure and possibly complete 
mission loss. Aerospace has been on the forefront in the development of high-fidelity analysis of 
composites, including the specialized tools that are required to evaluate delaminations and 
debonds in COPVs. 

Aerospace continues development of the COPV Stress Analysis Program (COSAP), begun in the 
1990s under Air Force and NASA sponsorship. COSAP is primarily comprised of pre-processing 
code that is designed to work with the commercial finite element code ABAQUS™. COSAP is 
flexible enough to model various vessel configurations, winding patterns and peculiar situations 
resulting from manufacturing considerations including delaminations. COSAP has been 
validated through experiments (Ref B- 17). Although COSAP is designed for analysis of COPVs, 
it may also be used for analysis of other filament-wound vessels, such as composite solid rocket 
motor cases. 

COSAP is designed to handle complicated features, such as the continuously changing local 
wrap angle and thickness resulting from the filament winding process. COSAP is designed to 
automatically generate an axisymmetric FEM for a COPV, while calculating the effective 3-D 
composite material properties. The input file requires major dimensions, basic material 
properties and definition of the winding pattern. The dome geometry can be specified by the user 
directly or optimized by the pre-processor. Depending upon the winding method, geodesic or 
planar, COSAP will automatically determine the local wrap angle and thicknesses of composite 
layers based on the input data and then generate the finite element mesh. To limit the total 
number of elements, a group of composite layers, sublaminate, can be incorporated into a row of 
elements. 

COSAP generates a FEM suitable for analysis with ABAQUS™. ABAQUS™ is used to perform 
a multi-step analysis that may include fabrication, sizing and burst cycles. Sizing, or autofrettage, 
is part of a typical manufacturing process, whereby the COPV is pressurized in order to yield the 
metal liner and as a result, the liner retains compressive residual stresses at lower operating 
pressure. This increases the number of cycles the tank can sustain before leaking. Upon 
completion of the multi-step analysis, the results are processed to evaluate burst pressure 
capability, fatigue life and the leak-before-burst failure mode. Prediction of the burst pressure is 
based upon progressive failure analysis described earlier, wherein for COPVs the ultimate failure 
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is usually fiber breakage. Preceding the global failure, delaminations may propagate into the boss 
region or into other regions of criticality, thereby leading to decreased fiber capability. 

The objective of this section is to demonstrate how COSAP or equivalent tools can be used to 
evaluate the burst strength reduction due to the presence of a delamination. A typical COPV 
with a 20-inch diameter and alternating hoop and helical plies with a total thickness of up to 
0.5 inches is analyzed, Figure B-40. In this example, the model development consisted of the 
following steps: 

The liner contour was found by digitizing the available drawings and using software to 
determine the exact shape; 

Ply thickness and initial ply angle at the tangent point were specified in the winding program, 
which then generated the remaining ply angle and thickness throughout; 

The modeling also considered the inclusion of proper ply drop off geometry and location; 
this data was obtained from a cross section of a constructed component; 

Candidate delamination locations were identified in critical ply interfaces, Figure B-41; 

With the above geometric information, COSAP was used to generate an axisymmetric FEM 
with cohesive elements located at candidate locations including the liner-to-overwrap 
interface. 



Figure B-40. A typical COPV with boundary conditions. 



Figure B-41 . Delamination location for this example. 


The analysis was performed in ABAQUS™ and included elastic-plastic material properties of 
the metallic liner and contact modeling was enforced to prevent delaminated surfaces from 
penetrating. Three load steps were required to properly evaluate the COPV: 1) loaded to 
autofrettage pressure of 4000 psi; 2) unloaded to 0 psi; and 3) loaded to maximum effective 
operating pressure, 3000 psi. Progressive failure analysis was active during these analysis steps. 
To build confidence in the model, strain gage data compared well to the predicted strains in the 
model. This was also true for the COPV axial and radial growth during autofrettage. 

Two failure modes that could result from an initial delamination were evaluated. The first failure 
mode is decreased burst capability due to a hypothetical fiber defect occurring directly above the 
delamination, Figure 42. The concern was that the delamination growth could amplify the effects 
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of local fiber damage. The fiber defect was modeled by substantially reducing the elastic 
modulus in the hoop direction. The analysis showed that while the fiber defect itself could 
substantially reduce the burst capability, the addition of a delamination near this hypothetical 
fiber damage reduced the strength by less than 5 percent. Based on previous experience the 
effects of delamination can be much larger than in this example, arising from multiple factors. 



Figure B-42. Hypothetical fiber defect near delamination. 

The second failure mode is delamination propagation towards the boss, Figure B-43. A 
delamination in this region could have uncertain effects due to the large external loads that 
exist in this region. Progressive failure analysis, as described earlier, was used to predict 
delamination growth. Fracture toughness properties were not available from test data, so values 
from 1. 0-5.0 lbs/in were assumed. The analysis predicted delamination growth in the region of 
ply drop-offs, which can be explained by the high transverse shear stresses present in this region. 
Once the delamination grew outside of the drop-off region, the delamination growth rate slowed 
substantially. Further, results showed that this growth was insensitive to the fracture toughness 
chosen, Figure B-44. Again, based on previous experience, delamination growth characteristics 
are driven by multiple factors and this example should not be used for the assessment of any 
other structure, especially considering that delamination growth to the boss is plausible. 



Figure B-43. Direction of delamination growth in this example. 
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Figure B-44. Predicted delamination growth for low and high fracture toughness values. 

In summary, an example was used to demonstrate how progressive failure analysis in 
conjunction with state-of-the-art pre-processing tools can be used to evaluate delaminations in 
COPVs. A limitation of this example is that not all failure modes were considered. While these 
tools can be used as part of an assessment of other failure modes, additional sophisticated 
analysis may be required. Particular to liner buckling, analysis tools need to be improved and 
validated through a test, non-destructive evaluation and analysis program. 
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